Do Correctional Facilities Correct Our Youth?: Effects of Incarceration and Court-Ordered Community Service on Personality Development
Kathryn L. Bollich-Ziegler
Emorie D. Beck
Patrick L. Hill
Joshua J. Jackson
15. January 2021
1 Workspace
1.1 Packages
1.2 Working Directory
1.3 Load Data
First, we will lead in a Google Sheet that indexes original variable names in the NLSY 1979 Child and Young Adults Sample as well as the original NLSY Sample of their mothers. The Sheet also contains rules splitting variables by type (matching variables, predictor variables, grouping variables, and outcome variables) as well as rules for compositing variables and their original scale.
# load codebook sheet containing variable names, recoding, etc #
codebook_path <- sprintf("%s/blob/main/data/codebook.xlsx?raw=true", data_path)
# create a meta file for downloading the child data
# dat %>% gs_read(ws = "Child") %>% select(`QName`, year) %>%
# write_delim(path = paste(data_path, "meta.CHILDYA", sep = "/"), delim = ",", col_names = F)
# create a meta file for downloading the mother data
# dat %>% gs_read(ws = "Mom") %>% select(`QName`, year) %>%
# write_delim(path = paste(data_path, "meta.NLSY79", sep = "/"), delim = ",", col_names = F)
p1f <- tempfile()
download.file(codebook_path, p1f, mode="wb")
# grab the child and mother data of variables chosen to include
child <- read_excel(p1f, sheet = "Child") %>% filter(Include == "Yes")
mom <- read_excel(p1f, sheet = "Mom") %>% filter(Include == "Yes")
match <- read_excel(p1f, sheet = "Items") %>% filter(Included == "x")
moder <- read_excel(p1f, sheet = "Moderators") # load in the Rdata files of the raw, uncleaned NLSY data sets
load(paste(local_path, "/data/nlsy79cya/nlsy79cya.RData", sep = "/"))
load(paste(local_path, "data/nlsy79/nlsy79moms.RData", sep = "/"))
# cut additional variables not included in the present analyses
new_data <- tbl_df(new_data) %>% select(one_of(child$RNUM))
mom_data <- tbl_df(mom_data) %>% select(one_of(mom$RNUM))
# get IDs of the moms to match with children
moms <- unique((new_data %>%
select(C0000100, C0000200) %>%
setNames(c("PROC_CID", "PROC_MID")))$PROC_CID)
# get IDs of the children
part_ids <- new_data %>%
select(C0000100, C0000200, C0005700, C0005400) %>%
setNames(c("PROC_CID", "PROC_MID", "Dem_DOBYear", "DemCSex")) %>%
filter(!is.na(PROC_CID) & !is.na(PROC_MID))
# find moms with children
moms <- unique(part_ids$PROC_MID)2 Data Cleaning
2.1 Matching Variables
Matching variables are those items that we will use to match individuals on. Ultimately, we will end up with a wide format data frame with the following variables:
knitr::kable(match %>% select(Category, Item, ItemText), "html", #"latex",
longtable = T, booktabs = T) %>%
# kable_styling(latex_options = c("striped","repeat_header"),full_width = F) %>%
kable_styling(bootstrap_options = c("striped"),full_width = F) %>%
#kable_styling(full_width = F) %>%
column_spec(1, width = "4cm") %>%
column_spec(3, width = "8cm") %>%
scroll_box(width = "750px", height = "400px")| Category | Item | ItemText |
|---|---|---|
| Background Variables | DemBirthOrder | Birth Order |
| Background Variables | Dem_DOBYear | Date of Birth |
| Background Variables | DemRace | Race |
| Background Variables | DemCSex | Sex |
| Birth & 1st Year of Life | BrnWkWnDueDate | Born in within 1 week of due date |
| Birth & 1st Year of Life | CSecBrth | Born via C-section |
| Birth & 1st Year of Life | DaysHosBrth | Days baby spent in hospital after birth |
| Birth & 1st Year of Life | Dr4Ill1stYr | Did child see a doctor in the 1st year of life for an illness |
| Birth & 1st Year of Life | BrthLnth | Length at birth |
| Birth & 1st Year of Life | BreastFed | Was Child Breastfed |
| Birth & 1st Year of Life | BrthWght | Weight at birth |
| Discipline | DscplnFreq | Composite of frequency of discipline |
| Discipline | Dscpln | Composite of groundings, lost allowance, being sent to rule, and lost privileges |
| Discipline | DSSpank | Use of corporal punishment/spanking |
| Early Childhood Health | HlthCndLim | Average Sum of Limiting Conditions in Childhood |
| Early Childhood Health | HlthCndLimAct | Child has a condition that limits activities |
| Early Childhood Health | HlthCndLimScl | Child has condition that limits school attendance |
| Early Childhood Health | HlthCndReqTr12mo | Child has health condition that required medical treatment in last 12 months |
| Early Childhood Health | HlthCndReqDr | Child has health condition that requires Dr. / treatment |
| Early Childhood Health | DemPheight | Height at age ? |
| Early Childhood Health | HlthMthrRtng | Mother’s rating of child health |
| Early Childhood Health | DemPweight | Weight at age ? |
| Early Delinquent / Antisocial Behavior | BPIAntiSoc | Behavioral Problems Index - Antisocial |
| Early Delinquent / Antisocial Behavior | BPIAnxDep | Behavioral Problems Index - Anxiety / Depression |
| Early Delinquent / Antisocial Behavior | BPIDepend | Behavioral Problems Index - Dependent |
| Early Delinquent / Antisocial Behavior | BPIHeadstrong | Behavioral Problems Index - Headstrong |
| Early Delinquent / Antisocial Behavior | BPIHyper | Behavioral Problems Index - Hyperactivity |
| Early Delinquent / Antisocial Behavior | BPIWdrawn | Behavioral Problems Index - Withdrawal |
| Early Delinquent / Antisocial Behavior | DelinqBeh | Composite of lying, staying out all night, staying out late, and skipping school |
| Early Delinquent / Antisocial Behavior | DelinqProp | Composite of stealing or vandalizing in the last year |
| Early Delinquent / Antisocial Behavior | DelinqDrnk | Drinking in the last year |
| Early Delinquent / Antisocial Behavior | DelinqVio | Harming someone bad enough to need doctor |
| Household Characteristics | HomeBin | Composite of Home binary valuables, including the presence or absence of books, TV, etc |
| Household Characteristics | ACT | Family activities, including eating together, reading together, going to museums, musicals, other outings |
| Household Characteristics | Home | Home “amenities”, including number of books, newspapers, musical instruments, music players, magazines |
| Household Characteristics | DemFathNHH | Type of Father in household – Bio Dad, Step Dad, No Dad |
| Mother Characteristics | FathGrade | Highest Grade of Father’s Education |
| Mother Characteristics | MEdHG | Highest Grade of Mom’s Education |
| Mother Characteristics | DemMAgeBirth | Mother’s age at birth |
| Mother Characteristics | MJail | Was mother ever in jail? |
| Pregnancy | MDiuretic | Did mother use a diuretic during pregnancy, due to doctor advice or otherwise |
| Pregnancy | MAlcohol | Did mother use alcohol during pregnancy |
| Pregnancy | MCoc | Did mother use cocaine during pregnancy |
| Pregnancy | MMarij | Did mother use Marijuana during pregnancy |
| Pregnancy | MAlcoholFreq | Frequency of maternal alcohol use during pregnancy |
| Pregnancy | MCocB4Preg | Maternal Cocaine Use before pregnancy |
| Pregnancy | MDrugLife | Maternal Drug Use at all |
| Pregnancy | MMarijB4Preg | Maternal Marijuana Use before pregnancy |
| Pregnancy | DemMVitamin | Mother Multivitamin during pregnancy |
| Psychological / Cognitive Characteristics & Attitudes | CogDgtBkwd | Cognitive Performance: Digits Backward |
| Psychological / Cognitive Characteristics & Attitudes | CogDgtFwd | Cognitive Performance: Digits Forward |
| Psychological / Cognitive Characteristics & Attitudes | CogPiatMath | Cognitive Performance: Piat Math |
| Psychological / Cognitive Characteristics & Attitudes | CogPiatRdCmp | Cognitive Performance: Piat Reading Reading Comprehension |
| Psychological / Cognitive Characteristics & Attitudes | CogPiatRdRec | Cognitive Performance: Piat Reading Reading Recognition |
| Psychological / Cognitive Characteristics & Attitudes | GenRoles | Gender Roles |
| Psychological / Cognitive Characteristics & Attitudes | EduRemed | Was child ever placed in remedial courses |
| Household Characteristics | HHJail | Was anyone in household ever in jail? |
| Early Childhood Health | HlthCndLimChck | How many of the checkbox health issues did the child have? |
2.1.1 Mom Data
# get vector of original reference numbers (column names) of the mom data
mom_long_cols <- (mom %>% filter(!is.na(Rule.Age)))$RNUM
# get vector of the new variable names of the mom data
mom_new_cols <- (mom %>% filter(!is.na(Rule.Age)))$comb_varname
mom_long <- mom_data %>%
filter(R0000100 %in% moms) %>%
select(R0000100, one_of(mom_long_cols)) %>% # select ID, DOB, and items for compositing
gather(key = item, value = value, one_of(mom_long_cols)) %>% # switch to long format
setNames(c("PROC_MID", "item", "value")) %>% # rename columns
full_join(part_ids) %>% # add participant IDs
mutate(item = mapvalues(item, mom_long_cols, mom_new_cols),
value = ifelse(grepl("MJail_MRes", item) == T & value == 5, 1,
ifelse(grepl("MJail_MRes", item) == T & !is.na(value) == T, 0,
ifelse(grepl("MJail_MRes", item) == F, value, NA))),
value = ifelse(grepl("ParDvrcd", item) == T & value == 3, 1,
ifelse(grepl("ParDvrcd", item) == T & !is.na(value) == T, 0,
ifelse(grepl("ParDvrcd", item) == F, value, NA))),
value = ifelse(grepl("HHMomDead", item) == T & value == 65, 1,
ifelse(grepl("HHMomDead", item) == T & !is.na(value) == T, 0,
ifelse(grepl("HHMomDead", item) == F, value, NA)))) %>% # recode items to to new names
left_join((mom %>% # add in rules for compositing and reverse coding
select(comb_varname, Rule.Age, Rule.Item, reverse_code) %>%
filter(comb_varname %in% mom_new_cols) %>%
rename(item = comb_varname))) 2.1.2 Child / YA Data
# create a small function for calculating the mode
Mode <- function(x) {
ux <- unique(x)
ux <- ux[!is.na(ux)]
ux[which.max(tabulate(match(x, ux)))]
}
# get vector of original reference numbers (column names) of the child data with multiple waves
child_long_cols <- (child %>% filter(year != "XRND" & grepl("match", var.type)))$RNUM
# get vector of the new variable names of the child data with multiple waves
child_new_cols <- (child %>% filter(year != "XRND" & grepl("match", var.type)))$comb_varname
match_long <- new_data %>%
# select ID, DOB, and items for compositing
select(C0000100, C0000200, C0005700, C0005400, one_of(child_long_cols)) %>%
gather(key = item, value = value, one_of(child_long_cols)) %>% # switch to long format
setNames(c("PROC_CID", "PROC_MID", "Dem_DOBYear", "DemCSex", "item", "value")) %>% # rename columns
mutate(item = mapvalues(item, child_long_cols, child_new_cols)) %>% # recode items to to new names
left_join((child %>%
select(comb_varname, Rule.Age, Rule.Item, reverse_code) %>% # add in rules for
filter(comb_varname %in% child_new_cols) %>% # compositing and reverse coding
rename(item = comb_varname))) %>%
full_join((mom_long)) %>%
separate(item, c("item", "year"), sep = "[.]") %>% # split items in item and year columns
mutate(year = as.numeric(year), # convert year to numeric
age = ifelse(!is.na(value) & year >= Dem_DOBYear,
year - Dem_DOBYear, NA)) # calculate age
match_long2 <- match_long %>%
full_join(match_long %>%
filter(grepl("HlthCndLimChck", item)) %>%
select(-item, -Rule.Age, -Rule.Item, -reverse_code) %>%
rename(hlthChck = value)) %>%
mutate(value = ifelse(grepl("HlthCndLim_", item) == T & is.na(value) == T, hlthChck, value)) %>%
# get rid of cases where we don't have an age or don't have a value / response
filter((!is.na(age) | !is.na(value)) & (age <= 14 |
(item %in% c("HHJail_HHJail", "HHMomDead_HHMomDead", "HHDadDead_HHDadDead")))) %>%
mutate(value = ifelse(grepl("UsualRes", item) == T,
mapvalues(value, seq(1,12,1), c(1, rep(0, 11))), # 1 = with mother; 2 = otherwise
ifelse(grepl("Limit", item) == T,
mapvalues(value, seq(3,0,-1), seq(0,3,1)), # reverse weird scale for Limiting / Structure
ifelse(grepl("HrTV", item) == T, value/100,
ifelse(grepl("HlthCndLim", item) == T & !is.na(value) == T & value != 0, 1,
ifelse(grepl("HlthCndLim", item) == T & !is.na(value) == T & value == 0, 0,
ifelse(grepl("HlthCndLimScl", item) == T & year %in% seq(1986, 2004,2),
mapvalues(value, seq(0,2,1), c(0,1,0)),
ifelse(grepl("EduRemed", item) == T & year %in% seq(1988,1998,2),
mapvalues(value, c(0,1,2,4), c(0,1,0,0)),
ifelse(item == "HlthCndReqTr12mo_HlthCndReqTr12mo",
mapvalues(value, seq(0,2,1), c(0,1,NA)),
ifelse(item == "DemFathNHH_FathNHH",
mapvalues(value, 1:4, c(1:3,1)),
ifelse(grepl("Mom2Dad", item) == T,
mapvalues(value, seq(1,13,1), c(1, rep(0,12))),
ifelse(grepl("FathGrade", item) == T,
mapvalues(value, c(seq(0,20,1),95), c(seq(0,20,1),NA)),
ifelse(grepl("Home", item) == T & value > 99, 99,
ifelse(grepl("NghdSf2Scl", item) == T,
mapvalues(value, c(seq(1,4,1), 6, 25), c(seq(1,4,1), NA, NA)),
value)))))))))))))) %>% #Recode weird time variable where 100 = 1 hr to 1 hr = 1 hr
group_by(item) %>%
mutate(min = min(value, na.rm = T), max = max(value, na.rm = T), # find min & max for reverse scoring
value = ifelse(reverse_code != "rev", value, # leave pos keyed items alone
reverse.code(-1, value, mini = min[1], maxi = max[1]))) %>% # reverse code neg keyed items
select(-min, -max) # get rid of min and max columns
seq_years <- seq(1988,1992, 2)
# create a sum variable for mom's response to tantrums #
match_long3 <- match_long2 %>%
filter(grepl("DSTantrum", item) == T) %>%
mutate(value = ifelse(grepl("DSTantrum", item) == T,
ifelse(year %in% seq_years & is.na(value) == F, 1, # different selected events in 88-92 coded as 1-7
ifelse(year %in% seq_years & is.na(value) == T, 0, # unselected events not coded
ifelse(year == 1996, recode(value, `2` = 0),
ifelse(grepl("Talk", item) == T, recode(value, `1` = 0, `0` = 1),
value)))), value)) %>%
group_by(PROC_CID, PROC_MID, year, age, item, Rule.Age, Rule.Item) %>%
summarize(value = max(value, na.rm = T), # some mothers filled out 6-9 and 10+ surveys in one year
value = ifelse(is.infinite(value) == T, NA, value)) %>%
group_by(PROC_CID, PROC_MID, year, age, Rule.Age, Rule.Item) %>% # compositing across different responses
summarize(value = sum(value, na.rm = T),
item = "MRTantrum_Tantrum") %>%
ungroup() %>%
# add in health data
full_join(match_long2 %>% filter(grepl("HlthCndLim_", item)) %>%
group_by(PROC_CID, PROC_MID, age, year, item, Rule.Age, Rule.Item) %>%
# average across unique items for the same age first
summarize(value = mean(value)) %>%
separate(item, c("type", "item")) %>%
group_by(PROC_CID, PROC_MID, age, year, type, Rule.Age, Rule.Item) %>%
# sum across broader categories of items for the same ages second
summarize(value = sum(value)) %>%
ungroup() %>%
mutate(item = paste(type, type, sep = "_"))) %>%
# add in all other match variables but the health and tantrum variables
full_join(match_long2 %>% filter(grepl("DSTantrum", item) == F & grepl("HlthCndLim_", item) == F))
# create a function for compositing matching variables
sum_fun <- function(df, Rule1, Rule2){
fun_call <- function(x, rule){
switch(rule,
average = mean(x, na.rm = T),
mode = Mode(x)[1],
sum = sum(x, na.rm = T),
select = unique(x)[1],
max = max(x, na.rm = T))
}
df %>%
group_by(PROC_CID, item, Dem_DOBYear) %>% # group by person and item (collapse across age)
summarize(value = fun_call(value, Rule1)) %>% # use proper summary function on each group combination
mutate(value = ifelse(is.nan(value) == T, NA, # get rid of NaN and Inf values
ifelse(is.infinite(value) == T, NA, value))) %>%
separate(item, c("type", "item")) %>% # split item further into category and unique item
group_by(PROC_CID, type) %>% # group by subject and type (collapse across item)
summarize(value = fun_call(value, Rule2)) %>% # use proper summary function on each group combination
mutate(value = ifelse(is.nan(value) == T, NA, # get rid of NaN and Inf values
ifelse(is.infinite(value) == T, NA, value))) %>%
ungroup()
}
child_wide_cols <- (child %>% filter(year == "XRND" & grepl("match", var.type)))$RNUM
child_new_cols <- (child %>% filter(year == "XRND" & grepl("match", var.type)))$comb_varname
child_new_cols <- gsub(".XRND", "", child_new_cols)
# here, we are going to grab the child data with only one response, clean up variables that need to be recoded
# these are intentionally in wide format and will eventually be joined with other data
match_wide <- new_data %>%
select(C0000100, C0000200, C0005700, C0005400, one_of(child_wide_cols)) %>% # select ID, DOB, and items for compositing
gather(key = item, value = value, -C0000100, -C0000200, -C0005700, -C0005400) %>% # switch to long format
setNames(c("PROC_CID", "PROC_MID", "Dem_DOBYear", "DemCSex", "item", "value")) %>% # rename columns
mutate(item = mapvalues(item, child_wide_cols, child_new_cols)) %>%
filter(!is.na(value) & !(item %in% c("DemHeightIN", "DemHeightFT", "DemWeightLBS"))) %>%
mutate(value = ifelse(grepl("MRed", item) == T, mapvalues(value, c(0,1,4), c(0,1,1)), value)) %>%
separate(item, c("type", "item"), sep = "[_]") %>% # split items in item and year columns
group_by(PROC_CID, PROC_MID, type, Dem_DOBYear, DemCSex) %>%
summarize(value = mean(value, na.rm = T)) %>%
ungroup() %>%
mutate(value = ifelse(grepl("MRed", type) == T | grepl("Diuret", type) == T, ceiling(value), value)) %>%
spread(key = type, value = value) %>%
mutate(DemRace = mapvalues(DemRace, seq(1,3,1), seq(3,1,-1)),
BreastFed = ifelse(WkBreastFedEnd > WkBreastFedBeg, WkBreastFedEnd - WkBreastFedBeg, NA)) %>%
select(-WkBreastFedEnd, -WkBreastFedBeg)
# run the summary function created above on the long format cross-time matching variables
match_wide <- match_long3 %>%
group_by(Rule.Age, Rule.Item) %>% # group by rules because it's faster
nest() %>% # nest all the data for each rule combination together
mutate(value = pmap(list(data, Rule.Age, Rule.Item), possibly(sum_fun, NA_real_))) %>% # call the summary function
select(-data) %>%
unnest(value) %>% # unnest the results
ungroup() %>%
select(-contains("Rule")) %>% # get rid of the rule variables
spread(key = type, value = value) %>% # change to wide format
mutate(HomeBin = ceiling(HomeBin),
MWelfare = ifelse(is.na(MWelfare) == T, 0, MWelfare),
HlthCndLim = ifelse(is.na(HlthCndLim) == T, 0, HlthCndLim)) %>%
select(-contains("eight")) %>% # get rid of height and weight; we'll take care of those later
full_join(match_wide) %>% # join with wide format matching variables
mutate_all(as.numeric) # make sure they are all numeric
## norm height and weight information by age) ##
# height and weight are tricky, so instead of using sometimes unreliable raw height and weight data, we will use
# the World Health Organization's norms to calculate height-for-age and weight-for-age percentiles
source(paste(data_path, "raw/main/data/who2007_R/who2007.r", sep = "/"))
htwt <- match_long3 %>% ungroup() %>%
# grab height and weight variables for 9-10 year olds
filter(item %in% c("DemHeightFT_HeightFT", "DemHeightIN_HeightIN", "DemWeightLBS_WeightLBS") & age %in% c(9,10)) %>%
select(PROC_CID, PROC_MID, age, DemCSex, year, item, value) %>%
group_by(PROC_CID, PROC_MID, age, DemCSex, year, item) %>%
summarize(value = mean(value, na.rm = T)) %>%
ungroup() %>%
spread(key = item, value = value) %>% # change to wide format
# get rid of implausible height values
filter(!(DemHeightFT_HeightFT > 5 & DemHeightIN_HeightIN < 12) &
!is.na(DemHeightIN_HeightIN) & !is.na(DemWeightLBS_WeightLBS) & !is.na(DemHeightFT_HeightFT)) %>%
# convert weight to kilograms, height to centimeters, and age to months
mutate(DemWeightKG_WeightKG = DemWeightLBS_WeightLBS/2.2,
DemHeightIN_HeightIN = ifelse(DemHeightIN_HeightIN < 12 & DemHeightFT_HeightFT < 6,
DemHeightIN_HeightIN + (DemHeightFT_HeightFT * 12), DemHeightIN_HeightIN),
DemHeightCM_HeightCM = DemHeightIN_HeightIN * 2.54,
ageMon = 12*age) %>%
# catch remaining implausible values
filter(DemWeightLBS_WeightLBS > 30 & DemHeightIN_HeightIN > 30)
htwt <- data.frame(htwt)
wfawho2007<-read.table(paste(data_path, "raw/main/data/who2007_R/wfawho2007.txt", sep = "/"), header=T,sep="",skip=0)
hfawho2007<-read.table(paste(data_path, "raw/main/data/who2007_R/hfawho2007.txt", sep = "/"),header=T,sep="",skip=0)
bfawho2007<-read.table(paste(data_path, "raw/main/data/who2007_R/bfawho2007.txt", sep = "/"),header=T,sep="",skip=0)
who2007(FilePath=paste(data_path, "data", sep = "/"),
FileLab = "nlsy_who", mydf=htwt, sex=DemCSex, age=ageMon,
weight=DemWeightKG_WeightKG, height=DemHeightCM_HeightCM)
htwt <- read.csv(paste(data_path, "raw/main/data/nlsy_who_z.csv", sep = "/"), stringsAsFactors = F) %>%
select(PROC_CID, PROC_MID, age, DemCSex, zhfa, zwfa) %>%
mutate(DemPweight = pnorm(zwfa), DemPheight = pnorm(zhfa)) %>%
select(-zwfa, -zhfa, -age) %>% ungroup()
# add the height and weight variables to the matching dataset
match_wide <- match_wide %>% left_join(htwt) %>% ungroup()
# get unique subjects in the matching variables
rm(list = ls()[grepl("match_long", ls())])2.2 Groups
To do propensity score weighting, we need to define our groups - Jail, Court Ordred Community Service, or neither. We are interested in incarceration in adolscence, which we define as ages below 19. Because we do not want our groups to overlap, anyone who was both incarcerated and participated in community service were included in the incarcerated group only.
# get vector of original reference numbers (column names) of the child data with multiple waves
child_long_cols <- (child %>% filter(year != "XRND" & var.type == "group"))$RNUM
# get vector of the new variable names of the child data with multiple waves
child_new_cols <- (child %>% filter(year != "XRND" & var.type == "group"))$comb_varname
# create age categories
age_orig <- seq(8, 19, 1)
age_cat <- rep(c(9,10,12,14,16,18), each = 2)
group_long <- new_data %>%
select(C0000100, C0000200, C0005700, C0005400, one_of(child_long_cols)) %>%
gather(key = item, value = value, one_of(child_long_cols)) %>%
setNames(c("PROC_CID", "PROC_MID", "Dem_DOBYear", "DemCSex", "item", "value")) %>% # rename columns
mutate(item = mapvalues(item, child_long_cols, child_new_cols)) %>%
separate(item, c("item", "year"), sep = "[.]") %>%
mutate(year = as.numeric(year),
age = ifelse(!is.na(value) & year >= Dem_DOBYear, year - Dem_DOBYear, NA)) %>%
filter(!is.na(age) & !is.na(value) & age <= 19) %>%
mutate(age = mapvalues(age, age_orig, age_cat),
#age = ifelse(age <= 19, "<19", ">19"),
value = ifelse(item == "Jail_ResJail",
mapvalues(value, seq(1:22), c(rep(0,4), 1, rep(0,17))), value)) %>%
#unite(item, item, age, sep = ".") %>%
group_by(PROC_CID, PROC_MID, DemCSex, item, age, year) %>%
summarize(value = max(value, na.rm = T)) %>%
ungroup() %>%
separate(item, c("type", "item")) %>%
group_by(PROC_CID, PROC_MID, DemCSex, type, year, age) %>%
summarize(value = max(value),
value = ifelse(value != 0, 1,
ifelse(is.na(value) == T, NA, 0)),
cont_age = ifelse(any(value == 1), min(age[value == 1], na.rm = T), NA),
cont_year = ifelse(any(value == 1), min(year[value == 1], na.rm = T), NA)) %>%
ungroup()
group_long2 <- group_long %>%
group_by(PROC_CID, PROC_MID, type) %>%
mutate_at(vars(cont_age, cont_year),
function(x){if(any(!is.na(x))) max(x, na.rm = T) else x}) %>%
group_by(PROC_CID, PROC_MID, DemCSex, type, cont_age, cont_year) %>%
summarize(value = max(value),
value = ifelse(value != 0, 1,
ifelse(is.na(value) == T, NA, 0))) %>%
ungroup()
group_wide <- group_long2 %>%
select(-cont_age, -cont_year) %>%
distinct() %>%
spread(key = type, value = value) %>%
gather(key = type, value = value, Cnsling:Jail) %>% # adding in NA's
mutate(value = mapvalues(value, c(0,1,NA), c(0,1,0))) %>%
spread(key = type, value = value) %>%
mutate(CommServ = ifelse(CrtCommServ == 1, 1, 0),
groups = ifelse(Jail == 1, 1,
ifelse(CommServ == 1, 2,
ifelse(Jail == 0 & CommServ == 0, 0, NA))),
groups = factor(groups, levels = seq(0,2,1), labels = c("None", "Jail", "CommServ")))
groups_tv <- group_long2 %>%
select(PROC_CID:DemCSex, groups = type, cont_age, cont_year) %>%
filter(complete.cases(.)) %>%
distinct() %>%
mutate(groups = ifelse(groups == "CrtCommServ", "CommServ", groups)) %>%
right_join(group_wide %>% select(-(Cnsling:Jail))) %>%
distinct() %>%
group_by(PROC_CID, PROC_MID, DemCSex, groups) %>%
summarize_at(vars(cont_age, cont_year), ~ifelse(all(is.na(.)), ., min(., na.rm = T))) %>%
ungroup()
match_wide <- match_wide %>%
left_join(group_wide %>% ungroup() %>% select(PROC_CID, groups))##
## None Jail CommServ
## 7098 399 2412.3 Personality
We have four predictor/outcomes, whose trajectory over time we will track using growth models. For these, we want data from adolescence, so we choose variables we choose 4 perosnality variables whose data we have throughout the period – Impulsivity, Sensation Seeking, Self-Esteem, and CESD depression.
child_long_cols <- (child %>% filter(year != "XRND" & var.type == "predictor"))$RNUM
child_new_cols <- (child %>% filter(year != "XRND" & var.type == "predictor"))$comb_varname
age_orig <- seq(6, 43, 1)
age_cat <- rep(seq(6,42,2), each = 2)
pred_alpha <- new_data %>%
select(C0000100, C0000200, C0005700, C0005400, one_of(child_long_cols)) %>%
gather(key = item, value = value, one_of(child_long_cols)) %>%
setNames(c("PROC_CID", "PROC_MID", "Dem_DOBYear", "DemCSex", "item", "value")) %>% # rename columns
mutate(item = mapvalues(item, child_long_cols, child_new_cols)) %>%
left_join((child %>%
select(comb_varname, reverse_code) %>% # add in rules for
filter(comb_varname %in% child_new_cols) %>% # compositing and reverse coding
rename(item = comb_varname))) %>%
separate(item, c("item", "year"), sep = "[.]") %>%
#filter(!grepl("chld", item)) %>%
mutate(year = as.numeric(year),
age = ifelse(!is.na(value) & year >= Dem_DOBYear, year - Dem_DOBYear, NA)) %>%
filter(!is.na(age) & !is.na(value) & age > 6) %>%
mutate(age = mapvalues(age, age_orig, age_cat)) %>%
#unite(item, item, age, sep = ".") %>%
group_by(item) %>%
mutate(value = ifelse(reverse_code != "rev", value,
reverse.code(-1, value, mini = 1, maxi = 4))) %>%
ungroup() %>%
mutate(item = gsub("chld", "", item)) %>%
group_by(PROC_CID, item, age, year) %>%
summarize(value = mean(value, na.rm = T)) %>%
ungroup()
pred_long <- pred_alpha %>%
separate(item, c("item", "type")) %>%
group_by(PROC_CID, item, age, year) %>%
summarize(value = mean(value, na.rm = T)) %>%
ungroup() %>%
mutate(age0 = age-14)
match_wide <- pred_long %>% filter(age == 14) %>%
select(PROC_CID, item, value) %>%
spread(item, value) %>%
right_join(match_wide)2.3.1 Scale Reliability
alpha_fun <- function(df){
(df %>% select(PROC_CID, type, value) %>%
spread(key = type, value = value) %>%
select(-PROC_CID) %>%
psych::alpha(.))$total[,1]
}
alpha.df <- pred_alpha %>%
separate(item, c("item", "type")) %>%
filter(type != "Capable") %>%
group_by(item, year) %>%
nest() %>%
mutate(alpha = map_dbl(data, alpha_fun))
alpha.tab <- alpha.df %>%
select(-data) %>%
spread(key = item, value = alpha) %>%
kable(., format = "html", booktabs = T, digits = 2,
caption = "<strong>Table S1</strong><br><em>Cronbach's Alpha Over Ages</br>") %>%
kable_styling(c("striped","repeat_header"),
full_width = F)
save_kable(alpha.tab, file = sprintf("%s/results/tables/Table_S%s_alphas.html", local_path, sup_tables))
sup_tables <<- sup_tables + 1
tables <- c(1)2.3.2 Raw Trajectories
First, we look at the raw trajectories of the data in order to ensure that our results match those from earlier work.
pred_long %>%
full_join(group_wide %>% select(PROC_CID, groups)) %>%
select(item, age, value, groups) %>%
filter(!is.na(item) & !is.na(groups)) %>%
Rmisc::summarySE(groupvars = c("item", "age", "groups"), measurevar = "value") %>%
group_by(item) %>%
#mutate(value = as.numeric(scale(value))) %>%
ggplot(aes(x = age, y = value)) +
geom_point(size = 1.3) +
geom_line(aes(color = groups)) +
#geom_errorbar(aes(ymin = value - ci, ymax = value + ci),position = "dodge", width = .1, stat = "identity") +
facet_wrap(~item, ncol = 2, scales = "free") +
theme_classic() +
theme(legend.position = "bottom")
2.4 Moderators
We test a variety of demographic moderators (e.g. gender) as well as moderators correpondeng to the ACE Model.
knitr::kable(moder %>% select(Category, Item, ItemText), "html", # "latex",
longtable = T, booktabs = T) %>%
# kable_styling(latex_options = c("striped","repeat_header"),full_width = F) %>%
kable_styling(bootstrap_options = c("striped","repeat_header"),full_width = F) %>%
#kable_styling(full_width = F) %>%
column_spec(1, width = "4cm") %>%
column_spec(3, width = "8cm")| Category | Item | ItemText |
|---|---|---|
| Household | HHMomDead | Is child’s mother deceased? |
| Household | HHDadDead | Is child’s father deceased? |
| Household | ParDvrcd | Are child’s parents divorced? |
| Household | DemFathNHH | Is child’s father in residence? |
| Demographic | DemCSex | Gender |
| Demographic | DemRace | Race |
| Physical Abuse | DSSpank | Was child spanked? |
| Emotional Abuse | ParEmoSupp | Emotional Support Raw Score |
| Health | HlthCndLim | Did child have any limiting health conditions? |
| SES | MWelfare | Was family ever on welfare? |
| Academic | SclSat | Satisfaction with School (Composite) |
| Academic | SchoolSat | Global satisfaction with school |
| Academic | SclChal | School is challenging (Composite) |
| Academic | EduRemed | Child was in remedial Math or English. |
moder_wide <- match_wide %>%
select(PROC_CID, Dem_DOBYear, one_of(moder$Item)) %>%
mutate(MWelfare = factor(ifelse(MWelfare == 0, 0,
ifelse(!is.na(MWelfare) == T, 1, NA)),
levels = c(0,1)),#, labels = c("NoWelfare","Welfare")),
DSSpank = factor(ifelse(DSSpank == 0, 0,
ifelse(!is.na(DSSpank) == T, 1, NA)),
levels = c(0,1)),#, labels = c("NotSpanked","Spanked")),
HlthCndLim = factor(ifelse(HlthCndLim > 1, 1,
ifelse(!is.na(HlthCndLim) == T, 0, NA)),
levels = c(0,1)),##, labels = c("NoHlthProb","HlthProb")),
DemRace = factor(mapvalues(DemRace, from = 1:3, to = 0:2),
levels = c(0,1,2)),#, labels = c("White", "Black", "Other")),
DemCSex = factor(mapvalues(DemCSex, from = 1:2, to = 0:1),
levels = c(0,1)),#, labels = c("M", "F")),
DemFathNHH = factor(mapvalues(DemFathNHH, from = 1:3, to = c(0,0,1)),
levels = 0:2),#, labels = c("Bio", "Step", "None")),
ParDvrcd = factor(ParDvrcd, levels = c(0,1)),#, labels = c("Tog", "Div")),
HHParDead = factor(Matrix::rowSums(cbind(HHMomDead, HHDadDead), na.rm = T), levels = c(0,1,2)),
HHDadDead = factor(HHDadDead, levels = c(0,1)),#, labels = c("DA", "DD")),
HHMomDead = factor(HHMomDead, levels = c(0,1)),#, labels = c("MA", "MD")),
none = NA_character_)
moder_long <- moder_wide %>% mutate_if(is.factor, funs(as.numeric(as.character(.)))) %>%
mutate_if(is.character, funs(as.numeric)) %>%
gather(key = moditem, value = modvalue, -PROC_CID, -Dem_DOBYear) %>%
filter(!is.na(modvalue))2.5 Outcomes
child_long_cols <- (child %>% filter(year != "XRND" & var.type == "outcome"))$RNUM
child_new_cols <- (child %>% filter(year != "XRND" & var.type == "outcome"))$comb_varname
age_orig <- seq(6, 25, 1)
age_cat <- rep(seq(6,24,2), each = 2)
outcome_long <- new_data %>%
select(C0000100, C0000200, C0005700, one_of(child_long_cols)) %>%
gather(key = item, value = value, one_of(child_long_cols)) %>%
setNames(c("PROC_CID", "PROC_MID", "Dem_DOBYear", "item", "value")) %>%
mutate(item = mapvalues(item, child_long_cols, child_new_cols)) %>%
separate(item, c("item", "year"), sep = "[.]") %>%
mutate(year = as.numeric(year),
age = ifelse(!is.na(value) & year >= Dem_DOBYear, year - Dem_DOBYear, NA)) %>%
filter(!is.na(age) & !is.na(value) & age <= 25) %>%
mutate(age = mapvalues(age, age_orig, age_cat)) %>%
#unite(item, item, age, sep = ".") %>%
group_by(PROC_CID, item, age) %>%
summarize(value = mean(value, na.rm = T)) %>%
ungroup()
outcome_wide <- outcome_long %>%
separate(item, c("item", "type")) %>%
group_by(PROC_CID, item, age) %>%
summarize(value = mean(value, na.rm = T)) %>%
ungroup() %>%
full_join(outcome_long) %>%
unite(comb, item, age, sep = ".") %>%
spread(key = comb, value = value)
save(match_wide, outcome_wide, pred_long, group_wide, moder_wide,
moder_long, moder, pred_alpha, groups_tv,
file = paste(data_path, "data", "data.RData", sep = "/"))3 Multiple Imputation
We perform multiple imputation using the mi package. We use the variables we defined before in the match_wide data set. Before doing so, we first create a missing data frame and examine missingness patterns. Then we run the mi and save them to a list for use in propensity score weighting.
3.1 Pre-Work: Correlations
r <- match_wide %>%
select(-PROC_CID, -PROC_MID) %>%
mutate_all(funs(as.numeric)) %>%
cor(., use = "pairwise")
raw.cor <- data.frame(r) %>% mutate(var = colnames(.)) %>%
gather(key = var2, value = value, -var) %>%
mutate(value = ifelse(var == var2, 0, value)) %>%
ggplot(aes(x = var, y = var2, fill = value)) +
geom_raster() +
labs(x = NULL, y = NULL, title = "Figure S1. Heatmap of Zero-Order Correlations Among Matching Variables") +
scale_fill_gradient2(low = "blue", mid = "white", high = "red",
limit = c(-1,1), midpoint = 0) +
theme_classic() +
theme(axis.text.x = element_text(angle = 90))
raw.cor
3.2 MI Procedure
outcome_subs <- unique(as.character(outcome_wide$PROC_CID))
match_subs <- unique(as.character(match_wide$PROC_CID))
group_subs <- unique(as.character(group_wide$PROC_CID))
pred_subs <- unique(as.character(pred_long$PROC_CID))
# figure out which subjects have grouping, matching, predictor, and outcome data
subs <- match_subs[match_subs %in% group_subs]
subs <- subs[subs %in% pred_subs]
subs <- subs[subs %in% outcome_subs]
# keep only those people
match_wide <- match_wide %>% filter(PROC_CID %in% subs)
pred_long <- pred_long %>% filter(PROC_CID %in% subs)
outcome_wide <- outcome_wide %>% filter(PROC_CID %in% subs)
# get a list of all objects in teh environment and keep only those we need to MI and PSW
all_env <- ls()
all_env <- all_env[!(all_env %in%
c("pred_long", "outcome_wide", "match_wide", "pred_alpha", "groups_tv",
"dat", "child", "mom", "data_path", "moder_wide", "group_wide",
"local_path", "sup_figs", "sup_tables", "moder_long"))]
rm(list = all_env)
gc()
# MI doesn't like tibbles, so we need to unclass and reclass the data
match_wide <- data.frame(unclass(match_wide)) # mi doesn't like tibbles
nlsy.mdf <- missing_data.frame(data.frame(unclass(match_wide)))
image(nlsy.mdf)# get descriptives of the data
des <- describe(match_wide,fast = T)
nlsy.mdf <- change(data = nlsy.mdf, y = rownames(des)[des$range >= 3],
what = "type", to = "continuous")
nlsy.mdf <- change(data = nlsy.mdf, y = c("DemPweight", "DemPheight"),
what = "type", to = "continuous")
gc()
nlsy.mi <- mi(nlsy.mdf, n.chains = 10, n.iter = 20, parallel = T)
gc()
image(nlsy.mi)completed.nlsy <- mi::complete(nlsy.mi)
save(completed.nlsy, nlsy.mdf, nlsy.mi, file = paste(data_path, "results/mi_dat.RData", sep = "/"))
clean.mi.fun <- function(df, Chain){
df <- tbl_df(df) %>% select(-contains("missing_"))
}
nested.mi <-
tibble(chain = paste("chain", 1:10, sep = "."), data = completed.nlsy) %>%
mutate(imp.data = map2(data, chain, possibly(clean.mi.fun, NA_real_)))
nested.mi <- nested.mi %>% mutate(match_set = "selection") %>%
bind_rows(nested.mi %>% mutate(match_set = "socialization"))4 Propensity Score Matching
Then, we perform propensity score weighting using our imputed datasets using the twang package. We then add the weights to our matching data frame, along with our predictor variables. To test the effectiveness of the propensity score weighting procedure, we examine the average standardized effect size in the balance tables. minimal effect sizes are candidates for beng dropped from the propensity score weighting, and large effect sizes mean our weighting procedure wasn’t effective. We can also examine these using balance plots.
4.1 Raw Group Differences
temp <- tbl_df(match_wide) %>%
group_by(groups) %>%
summarize(n = n(),
`% Female` = sum(DemCSex == 2)/ sum(DemCSex %in% c(1,2)),
`Mean Age 2017` = mean(2017 - Dem_DOBYear, na.rm = T))
age_tab <- bind_rows(temp, temp %>% ungroup() %>%
summarize(n = sum(n),
`% Female` = mean(`% Female`), groups = "Total",
`Mean Age 2017` = mean(`Mean Age 2017`))) %>%
kable(., "html", #"latex",
escape = F, booktabs = T, digits = 2,
col.names = c("Group", "N", "% Female", "Mean Age 2017")
, caption = "<strong>Table S2</strong><br><em>Count, Age, and Gender Breakdown Across Groups</em>") %>%
# kable_styling(latex_options = c("striped","repeat_header"),full_width = F)
kable_styling(bootstrap_options = c("striped","repeat_header"),full_width = F)
save_kable(age_tab, file = sprintf("%s/results/tables/Table_S%s_descriptives.html", local_path, sup_tables))
sup_tables <<- sup_tables + 1
age_tab| Group | N | % Female | Mean Age 2017 |
|---|---|---|---|
| None | 7096 | 0.51 | 30.60 |
| Jail | 399 | 0.24 | 31.20 |
| CommServ | 241 | 0.41 | 29.44 |
| Total | 7736 | 0.39 | 30.41 |
4.2 Choose Matching Variables
We had some issues achieving balance. So, we elected to evaluate the effect size of our matching variables across groups and keep only those variables with raw differences larger than .2.
df <- tbl_df(completed.nlsy[[1]]) %>% select(-contains("missing"))
df <- df %>% select(-groups) %>%
mutate_if(is.factor, funs(as.numeric(as.character(.)))) %>%
full_join(df %>% select(PROC_CID, groups)) %>%
gather(key = item, value = value, CESD:DemPheight)
d_fun <- function(Item, V1, V2){
dat <- df %>% filter(item == Item & groups %in% c(V1, V2)) %>%
mutate(groups = factor(groups, levels = c(V1, V2)))
d <- lsr::cohensD(value ~ groups, data = dat)
}
keep <- crossing(
Groups = c("None-Jail", "None-CommServ", "Jail-CommServ"),
Item = unique(df$item)
) %>%
tbl_df() %>%
separate(Groups, c("V1", "V2")) %>%
mutate(d = pmap_dbl(list(Item, V1, V2), d_fun))
keep_vars <- unique((keep %>% filter(d > .2))$Item)
keep_vars <- c(keep_vars, "SelfEstm")4.3 Define Functions
# this function actually runs the propensity score weighting procedure
psw_fun <- function(df, match_set){
if(match_set == "socialization"){
to.match <- colnames(df)[-which(colnames(df) %in% c("PROC_CID","PROC_MID","groups"))]
} else{
to.match <- colnames(df)[-which(colnames(df) %in% c("PROC_CID","PROC_MID","groups",
"CESD", "SelfEstm", "Impls", "SensSeek"))]
}
match.formula <- as.formula(paste("groups ~ ", paste(to.match, collapse=" + "), sep = " "))
y <- mnps(formula = match.formula,
data=df,
n.trees=10000,
stop.method = c("es.mean"))
}
# changing the data fed into psw into a data frame because it won't work with tibbles
psw_df <- function(psw){psw$data <- data.frame(psw$data); psw}
# once weighting is done, we pull the weights back into the original mi df
# and join back with the df of personality variables
add_psw_weights <- function(psw, df){
df$psw.weights <- get.weights(psw)
df <- df %>%
select(PROC_CID, groups, psw.weights, CESD, SelfEstm, SensSeek, Impls) %>%
mutate(groups = relevel(groups, "None"))
return(df)
}
# this function creates the balance table of the psw weights and filters
# the results into variables the matching procedure did not fix and
# those that it did
unbalanced_fun <- function(x){
#x <- bal.table(psw)
y <- x %>%
filter(stop.method == "ks.mean" | stop.method == "unw") %>%
select(stop.method, tmt1, tmt2, var, std.eff.sz) %>%
group_by(stop.method, tmt1, tmt2, var) %>%
summarize(m.eff = mean(std.eff.sz)) %>%
filter((stop.method == "ks.mean" & m.eff >= .2) |
(stop.method == "unw" & m.eff <= .05)) %>%
setNames(c("stop.method", "Group1", "Group2", "Variable", "Average Standardized Effect"))
unbalanced.var <- tbl_df(y) %>% filter(stop.method == "ks.mean") %>% select(-stop.method)
smalldiff.var <- tbl_df(y) %>% filter(stop.method == "unw") %>% select(-stop.method)
return(list(unbalanced = unbalanced.var,smalldiff = smalldiff.var))
}4.4 Run Matching
plan(multiprocess(workers = 3L))
nested.psw <- nested.mi %>%
mutate(imp.data = map(imp.data, function(x){x %>% select(PROC_CID, groups, one_of(keep_vars))}),
psw = future_map2(imp.data, match_set, possibly(psw_fun, NA_real_), .progress = T),
psw = map(psw, possibly(psw_df, NA_real_)),
bal.tab = map(psw, possibly(bal.table, NA_real_)),
bal.df = map(bal.tab, possibly(unbalanced_fun, NA_real_)),
unbal.tab = map(bal.df, possibly(~.$unbalanced, NA_real_)),
smalldiff.tab = map(bal.df, possibly(~.$smalldiff, NA_real_)))
closeAllConnections()
gc()4.5 Balance Tables
Once propensity scores are estimated, bal.table() produces a table that shows how well the resulting weights succeed in manipulating the groups so that they match on pre-adolescent matching characteristics.
load(sprintf("%s/results/psw.RData", local_path))
# this table identifies variables differences in effect size pre and post weighting
diff.tab <- nested.psw %>% filter(match_set == "socialization") %>%
unnest(bal.tab) %>%
filter(stop.method %in% c("es.mean", "unw")) %>%
group_by(stop.method, tmt1, tmt2, var) %>%
summarise(m.eff = mean(std.eff.sz)) %>%
select(tmt1, tmt2, var, stop.method, m.eff) %>%
spread(key = stop.method, value = m.eff) %>%
mutate(diff = unw - es.mean) %>%
arrange(diff)
diff.tab <- kable(diff.tab, "html", digits = 3,
longtable = T, booktabs = T, escape = F,
col.names = c("Group 1", "Group 2", "Matching<br>Variable",
"Weighted<br>Effect Size", "Unweighted<br>Effect Size",
"Change<br>Effect Size"),
caption = "<strong>Table S3</strong><br><em>Change in Effect Size after Propensity Score Weighting<em>") %>%
# kable_styling(latex_options = c("striped","repeat_header"),full_width = F)
kable_styling(bootstrap_options = c("striped","repeat_header"),full_width = F) %>%
scroll_box(width = "750px", height = "800px")
save_kable(diff.tab, file = sprintf("%s/results/tables/Table_S%s_PSW_diff.html", local_path, sup_tables))
diff.tab| Group 1 | Group 2 |
Matching Variable |
Weighted Effect Size |
Unweighted Effect Size |
Change Effect Size |
|---|---|---|---|---|---|
| Jail | CommServ | RelDad | 0.101 | 0.031 | -0.070 |
| Jail | CommServ | MNetInc | 0.079 | 0.023 | -0.055 |
| Jail | CommServ | ParEmoSupp | 0.084 | 0.030 | -0.054 |
| Jail | CommServ | RelTmFam | 0.135 | 0.081 | -0.053 |
| Jail | CommServ | MDrugLife | 0.160 | 0.117 | -0.043 |
| Jail | CommServ | PeerPress | 0.097 | 0.063 | -0.034 |
| Jail | CommServ | MNbhd | 0.116 | 0.082 | -0.034 |
| Jail | CommServ | SelfEstm | 0.099 | 0.068 | -0.031 |
| Jail | CommServ | EduRemed:0 | 0.034 | 0.013 | -0.020 |
| Jail | CommServ | EduRemed:1 | 0.034 | 0.013 | -0.020 |
| Jail | CommServ | Delinq | 0.057 | 0.042 | -0.015 |
| None | Jail | MRTantrum | 0.073 | 0.059 | -0.014 |
| Jail | CommServ | SensSeek | 0.041 | 0.029 | -0.012 |
| None | CommServ | MJail:0 | 0.030 | 0.023 | -0.007 |
| None | CommServ | MJail:1 | 0.030 | 0.023 | -0.007 |
| None | Jail | MCocB4Preg:0 | 0.052 | 0.050 | -0.002 |
| None | Jail | MCocB4Preg:1 | 0.052 | 0.050 | -0.002 |
| None | Jail | Dem_DOBYear | 0.104 | 0.104 | -0.001 |
| None | CommServ | DemMAgeBirth | 0.158 | 0.162 | 0.004 |
| Jail | CommServ | HHJail:0 | 0.063 | 0.069 | 0.006 |
| Jail | CommServ | HHJail:1 | 0.063 | 0.069 | 0.006 |
| Jail | CommServ | MWelfare | 0.032 | 0.065 | 0.033 |
| None | CommServ | Dem_DOBYear | 0.167 | 0.201 | 0.034 |
| Jail | CommServ | SclSat | 0.051 | 0.085 | 0.034 |
| Jail | CommServ | MCocB4Preg:0 | 0.071 | 0.105 | 0.034 |
| Jail | CommServ | MCocB4Preg:1 | 0.071 | 0.105 | 0.034 |
| None | CommServ | ParDvrcd:0 | 0.105 | 0.144 | 0.038 |
| None | CommServ | ParDvrcd:1 | 0.105 | 0.144 | 0.038 |
| None | CommServ | GenRoles | 0.050 | 0.092 | 0.042 |
| None | Jail | RelDad | 0.163 | 0.206 | 0.043 |
| Jail | CommServ | Mom2Dad:0 | 0.183 | 0.228 | 0.044 |
| Jail | CommServ | Mom2Dad:1 | 0.183 | 0.228 | 0.044 |
| None | CommServ | MDrugLife | 0.174 | 0.220 | 0.046 |
| Jail | CommServ | MSmoke:0 | 0.097 | 0.144 | 0.046 |
| Jail | CommServ | MSmoke:1 | 0.097 | 0.144 | 0.046 |
| Jail | CommServ | Dscpln | 0.080 | 0.128 | 0.047 |
| Jail | CommServ | SclChal | 0.041 | 0.089 | 0.048 |
| Jail | CommServ | GenRoles | 0.040 | 0.092 | 0.052 |
| None | Jail | SensSeek | 0.152 | 0.210 | 0.057 |
| None | Jail | RelTmFam | 0.115 | 0.179 | 0.064 |
| Jail | CommServ | CESD | 0.039 | 0.105 | 0.066 |
| None | CommServ | SelfEstm | 0.074 | 0.147 | 0.073 |
| None | CommServ | MSmoke:0 | 0.029 | 0.104 | 0.075 |
| None | CommServ | MSmoke:1 | 0.029 | 0.104 | 0.075 |
| None | CommServ | MRTantrum | 0.118 | 0.200 | 0.082 |
| None | Jail | MDrugLife | 0.020 | 0.103 | 0.083 |
| Jail | CommServ | Impls | 0.076 | 0.162 | 0.086 |
| None | CommServ | MCocB4Preg:0 | 0.036 | 0.125 | 0.089 |
| None | CommServ | MCocB4Preg:1 | 0.036 | 0.125 | 0.089 |
| None | CommServ | PeerPress | 0.063 | 0.153 | 0.090 |
| Jail | CommServ | MRTantrum | 0.050 | 0.141 | 0.091 |
| None | Jail | SclChal | 0.013 | 0.108 | 0.094 |
| None | Jail | HHJail:0 | 0.126 | 0.221 | 0.095 |
| None | Jail | HHJail:1 | 0.126 | 0.221 | 0.095 |
| None | CommServ | SensSeek | 0.125 | 0.227 | 0.101 |
| None | Jail | DemBirthOrder | 0.050 | 0.156 | 0.106 |
| None | Jail | CESD | 0.075 | 0.185 | 0.110 |
| None | CommServ | Mom2Dad:0 | 0.081 | 0.193 | 0.112 |
| None | CommServ | Mom2Dad:1 | 0.081 | 0.193 | 0.112 |
| None | CommServ | MNbhd | 0.097 | 0.216 | 0.119 |
| None | CommServ | SclSat | 0.037 | 0.158 | 0.121 |
| None | CommServ | MNetInc | 0.057 | 0.178 | 0.122 |
| Jail | CommServ | BPI | 0.050 | 0.175 | 0.126 |
| Jail | CommServ | ParDvrcd:0 | 0.054 | 0.188 | 0.134 |
| Jail | CommServ | ParDvrcd:1 | 0.054 | 0.188 | 0.134 |
| None | Jail | MJail:0 | 0.010 | 0.149 | 0.138 |
| None | Jail | MJail:1 | 0.010 | 0.149 | 0.138 |
| Jail | CommServ | DemBirthOrder | 0.038 | 0.178 | 0.140 |
| None | CommServ | Impls | 0.032 | 0.184 | 0.152 |
| Jail | CommServ | MJail:0 | 0.020 | 0.173 | 0.153 |
| Jail | CommServ | MJail:1 | 0.020 | 0.173 | 0.153 |
| None | Jail | GenRoles | 0.029 | 0.185 | 0.156 |
| None | Jail | SelfEstm | 0.038 | 0.195 | 0.157 |
| None | Jail | DemMAgeBirth | 0.025 | 0.186 | 0.161 |
| None | CommServ | RelDad | 0.067 | 0.230 | 0.163 |
| None | CommServ | SclChal | 0.029 | 0.196 | 0.167 |
| None | CommServ | DemCSex:1 | 0.040 | 0.212 | 0.172 |
| None | CommServ | DemCSex:2 | 0.040 | 0.212 | 0.172 |
| None | Jail | MNetInc | 0.028 | 0.202 | 0.174 |
| None | Jail | SclSat | 0.066 | 0.243 | 0.177 |
| None | Jail | MSmoke:0 | 0.070 | 0.247 | 0.177 |
| None | Jail | MSmoke:1 | 0.070 | 0.247 | 0.177 |
| None | CommServ | EduRemed:0 | 0.028 | 0.209 | 0.181 |
| None | CommServ | EduRemed:1 | 0.028 | 0.209 | 0.181 |
| None | Jail | PeerPress | 0.035 | 0.215 | 0.181 |
| None | Jail | ParDvrcd:0 | 0.146 | 0.331 | 0.185 |
| None | Jail | ParDvrcd:1 | 0.146 | 0.331 | 0.185 |
| None | Jail | EduRemed:0 | 0.029 | 0.221 | 0.192 |
| None | Jail | EduRemed:1 | 0.029 | 0.221 | 0.192 |
| Jail | CommServ | DemMAgeBirth | 0.145 | 0.348 | 0.203 |
| None | CommServ | CESD | 0.072 | 0.288 | 0.216 |
| None | CommServ | HHJail:0 | 0.064 | 0.290 | 0.226 |
| None | CommServ | HHJail:1 | 0.064 | 0.290 | 0.226 |
| None | CommServ | BPI | 0.040 | 0.273 | 0.234 |
| Jail | CommServ | DemCSex:1 | 0.138 | 0.376 | 0.238 |
| Jail | CommServ | DemCSex:2 | 0.138 | 0.376 | 0.238 |
| None | CommServ | RelTmFam | 0.021 | 0.260 | 0.239 |
| Jail | CommServ | Dem_DOBYear | 0.063 | 0.305 | 0.242 |
| None | Jail | MNbhd | 0.026 | 0.299 | 0.272 |
| None | Jail | ParEmoSupp | 0.028 | 0.309 | 0.281 |
| None | CommServ | ParEmoSupp | 0.056 | 0.339 | 0.284 |
| None | CommServ | DemBirthOrder | 0.048 | 0.334 | 0.286 |
| None | Jail | Impls | 0.056 | 0.346 | 0.290 |
| None | Jail | MWelfare | 0.024 | 0.319 | 0.295 |
| None | Jail | Mom2Dad:0 | 0.113 | 0.421 | 0.308 |
| None | Jail | Mom2Dad:1 | 0.113 | 0.421 | 0.308 |
| None | Jail | BPI | 0.089 | 0.449 | 0.359 |
| None | CommServ | MWelfare | 0.016 | 0.384 | 0.368 |
| None | CommServ | Dscpln | 0.014 | 0.404 | 0.390 |
| None | Jail | DemCSex:1 | 0.178 | 0.587 | 0.409 |
| None | Jail | DemCSex:2 | 0.178 | 0.587 | 0.409 |
| None | Jail | Dscpln | 0.080 | 0.531 | 0.451 |
| None | Jail | Delinq | 0.110 | 0.572 | 0.462 |
| None | CommServ | Delinq | 0.053 | 0.610 | 0.557 |
sup_tables <<- sup_tables + 1
# this table shows variables that are not matched after weighting
unbal.tab <- nested.psw %>%
filter(match_set == "socialization") %>%
select(chain, match_set, unbal.tab) %>%
unnest(unbal.tab) %>%
group_by(Group1, Group2, Variable) %>%
summarize(`Average Standardized Effect` = mean(`Average Standardized Effect`, na.rm = T)) %>%
ungroup()
# kable(unbal.tab, "latex", longtable = T, booktabs = T,
# caption = "Unbalanced Variables after Propensity Score Weighting") %>%
# kable_styling(latex_options = c("striped","repeat_header"),full_width = F)
# this table shows variables that were already matched prior to weighting
smalldiff.tab <- nested.psw %>%
filter(match_set == "socialization") %>%
select(chain, match_set, smalldiff.tab) %>%
unnest(smalldiff.tab) %>%
group_by(Group1, Group2, Variable) %>%
summarize(`Average Standardized Effect` = mean(`Average Standardized Effect`, na.rm = T)) %>%
ungroup()
kable(smalldiff.tab
, "html" # "latex",
, longtable = T
, booktabs = T
, caption = "Variables with Initial Small Difference"
, digits = 2) %>%
# kable_styling(latex_options = c("striped","repeat_header"),full_width = F)
kable_styling(bootstrap_options = c("striped","repeat_header"),full_width = F) %>%
scroll_box(width = "750px", height = "400px")| Group1 | Group2 | Variable | Average Standardized Effect |
|---|---|---|---|
| Jail | CommServ | CESD | 0.02 |
| Jail | CommServ | Delinq | 0.03 |
| Jail | CommServ | Dscpln | 0.05 |
| Jail | CommServ | EduRemed:0 | 0.01 |
| Jail | CommServ | EduRemed:1 | 0.01 |
| Jail | CommServ | GenRoles | 0.04 |
| Jail | CommServ | HHJail:0 | 0.05 |
| Jail | CommServ | HHJail:1 | 0.05 |
| Jail | CommServ | MCocB4Preg:0 | 0.04 |
| Jail | CommServ | MCocB4Preg:1 | 0.04 |
| Jail | CommServ | MNetInc | 0.02 |
| Jail | CommServ | ParEmoSupp | 0.03 |
| Jail | CommServ | PeerPress | 0.03 |
| Jail | CommServ | RelDad | 0.03 |
| Jail | CommServ | SclChal | 0.04 |
| Jail | CommServ | SclSat | 0.02 |
| Jail | CommServ | SelfEstm | 0.02 |
| Jail | CommServ | SensSeek | 0.02 |
| None | CommServ | MCocB4Preg:0 | 0.00 |
| None | CommServ | MCocB4Preg:1 | 0.00 |
| None | CommServ | MJail:0 | 0.02 |
| None | CommServ | MJail:1 | 0.02 |
| None | Jail | MCocB4Preg:0 | 0.02 |
| None | Jail | MCocB4Preg:1 | 0.02 |
4.6 Balance Plots
In these plots, substantial reductions in effect sizes are observed for most variables (blue lines), with only one variable showing an increase in effect size (red lines), but only a seemingly trivial increase. Closed red circles indicate a statistically significant difference, many of which occur before weighting, none after.
5 Selection Effects: Logistic Regressions
First, we want to test selection effects – does personality predict whether an adolescent has contact with the criminal justice system. To that, we use personality at age 14 in a multinomial regression using the multinom() function in the nnet package.
The propensity score weights we use here are different from the ones we’ll use for socialization effects because they do not account for baseline personality when matching. ## Run Models
load(sprintf("%s/results/psw_small.RData", data_path))
nested.sel <- expand.grid(
Trait = unique(pred_long$item),
weight = c("weighted", "unweighted"),
stringsAsFactors = F
)
sel_fun <- function (trait, weight, ms){
models <- lapply(1:10, function(x){
# require(tidyverse); require(lavaan)
print(paste(trait, x), sep = " ")
df <- nested.psw %>% filter(chain == paste("chain", x, sep = ".") & match_set == ms) %>%
unnest(imp.data.long) %>%
left_join(pred_long %>% filter(age == 14 & item == trait)) %>%
filter(!is.na(item))
if(weight != "weighted"){
df$psw.weights <- rep(1, nrow(df))
}
mod <- nnet::multinom(
groups ~ value, data = df, weights = psw.weights
)
# file <- sprintf("%s/results/models/%s_%s_chain%s.RData", data_path, trait, event, x)
# save(mod, file = file)
return(mod)
})
}
start <- Sys.time()
nested.sel <- nested.sel %>%
mutate(sel.mod = map2(Trait, weight, ~sel_fun(.x, .y, "selection")))
end <- Sys.time()
save(nested.sel, file = sprintf("%s/results/selection.RData", data_path))5.1 Pool Results
load(url(sprintf("%s/blob/main/results/selection.RData?raw=true", data_path)))
sel_pool_fun <- function(models, type){
###
### Get coefs for pooling FE's across imputations ###
####
nchains <- length(models)
nterms <- nrow(vcov(models[[1]]))
names <- crossing(
group = rownames(summary(models[[1]])$coefficients),
term = colnames(summary(models[[1]])$coefficients)
) %>% arrange(term, rev(group)) %>%
unite(term, group, term, sep = ":")
coefs <- sapply(models, function(x) summary(x)$coefficients)
variances <- sapply(models, FUN = function(x) diag(vcov(x)))
rownames(coefs) <- names$term; rownames(variances) <- rownames(vcov(models[[1]]))
variances <- variances[names$term,]
# average effects for each term
coefs_mean <- apply(coefs, 1, mean)
var_mean <- apply(variances, 1, mean)
# variance of fixed effects
coefs_var <- apply(coefs, 1, var)
e <- t(coefs) - matrix(coefs_mean, nrow = nchains, ncol = nterms, byrow = TRUE)
b <- (t(e) %*% e)/(nchains - 1)
t <- var_mean + (1 + 1/nchains) * b
r <- (1 + 1/nchains) * diag(b/var_mean)
lambda <- (1 + 1/nchains) * diag(b/t)
dfcom <- df.residual(models[[1]])
if (is.null(dfcom)) {
dfcom <- 999999
warning("Large sample assumed.")
}
lambda[lambda < 1e-04] <- 1e-04
dfold <- (nchains - 1)/lambda^2
dfobs <- (dfcom + 1)/(dfcom + 3) * dfcom * (1 - lambda)
df <- dfold * dfobs/(dfold + dfobs)
se <- sqrt(diag(t))
t.val <- coefs_mean/se
p <- 2 * (1 - pt(abs(t.val), df = df))
CI = qt(.975, df = df)*se
fmi <- (r + 2/(df + 3))/(r + 1)
names(r) <- names(df) <- names(fmi) <- names(lambda) <- names$term
fixeff <- tibble(
term = names$term,
Estimate = coefs_mean,
# t = t.val,
CI.lower = coefs_mean-CI,
CI.upper = coefs_mean+CI
)
}5.2 Table Results
nested.sel <- nested.sel %>%
mutate(pool = map(sel.mod, possibly(sel_pool_fun, NA_real_)))
sel.tab <- nested.sel %>%
select(-sel.mod) %>%
unnest(pool, .drop = T) %>%
mutate(sig = ifelse(sign(CI.lower) == sign(CI.upper), "sig", "ns")) %>%
separate(term, c("group", "term"), sep = ":") %>%
mutate(group = mapvalues(group, c("Jail", "CommServ"), c("Incarcerated", "COV")),
Trait = mapvalues(Trait, traits$old, traits$breaks),
Trait = factor(Trait, levels = traits$breaks)) %>%
arrange(term) %>%
mutate_at(vars(Estimate, CI.lower, CI.upper), funs(exp)) %>%
filter(term == "value")
save(sel.tab, file = sprintf("%s/results/selection_res.RData", local_path))
sel.kable <- sel.tab %>%
mutate(est = sprintf("%.2f [%.2f, %.2f]", Estimate, CI.lower, CI.upper)
, est = ifelse(sig == "sig", sprintf("<strong>%s</strong>", est), est)) %>%
select(Trait, weight, group, est) %>%
pivot_wider(names_from = c("group", "weight")
, values_from = "est") %>%
select(Trait, Incarcerated_unweighted, Incarcerated_weighted, COV_unweighted, COV_weighted) %>%
arrange(Trait) %>%
kable(.
, "html"
, escape = F
, caption = "<strong>Table 1</strong><br><em>Unweighted and Weighted Effects of Baseline Personality on the Probability of Experiencing Incarceration or Community Service</em>"
, col.names = c("Trait", rep("OR [CI]", 4))
, align = c("r", rep("c", 4))) %>%
kable_styling(full_width = F) %>%
add_header_above(c(" " = 1, "Unweighted" = 1, "Weighted" = 1, "Unweighted" = 1, "Weighted" = 1)) %>%
add_header_above(c(" " = 1, "Incarceration" = 2, "Community Service" = 2)) %>%
add_footnote(" The reference group was individuals neither in the Incarceration nor Community Service groups. OR = Odds ratio. CI = Confidence interval. Bolded values represent estimated ORs whose CIs do not overlap with 1.")
sel.kable| Trait | OR [CI] | OR [CI] | OR [CI] | OR [CI] |
|---|---|---|---|---|
| Sensation Seeking | 1.39 [1.16, 1.67] | 1.54 [1.30, 1.84] | 1.43 [1.14, 1.80] | 1.41 [1.16, 1.71] |
| Impulsivity | 1.81 [1.49, 2.20] | 1.37 [1.19, 1.58] | 1.46 [1.15, 1.86] | 1.07 [0.87, 1.30] |
| Depressive Symptoms | 1.48 [1.13, 1.95] | 1.64 [1.39, 1.94] | 1.67 [1.20, 2.32] | 1.48 [1.09, 2.02] |
| Self- Esteem | 0.59 [0.41, 0.87] | 0.81 [0.56, 1.15] | 0.63 [0.39, 1.01] | 1.20 [0.82, 1.76] |
| a The reference group was individuals neither in the Incarceration nor Community Service groups. OR = Odds ratio. CI = Confidence interval. Bolded values represent estimated ORs whose CIs do not overlap with 1. |
save_kable(sel.kable, file = sprintf("%s/results/tables/Table_1_Selection.html", local_path))
sel.tab %>%
filter(weight == "weighted") %>%
ggplot(aes(x = group, y = Estimate)) +
# scale_y_continuous(limits = c(-1,2), breaks = seq(-3,3,1))+
scale_color_manual(values = c("springgreen3", "blue")) +
geom_hline(aes(yintercept = exp(0)), linetype = "dashed") +
geom_errorbar(aes(ymin = CI.lower, ymax = CI.upper), width = .2, size = .75) +
geom_point(aes(color = group), size = 4) +
labs(y = "Odds", x = NULL) +
# coord_flip() +
facet_grid(.~Trait, scale = "free") +
theme_classic() +
theme(legend.position = "none",
axis.text = element_text(face = "bold", size = rel(1.2)),
axis.text.x = element_text(face = "bold", size = rel(1.2), hjust = 1, angle = 45),
axis.title = element_text(face = "bold", size = rel(1.2)),
strip.text = element_text(face = "bold", size = rel(1.2)))
7 Moderators
7.1 Functions
7.1.1 Model Functions
mod_fun2 <- function(dv, Moder, Form, groups, psw.weights){
print(paste(dv, Moder, Form, sep = " "))
mods <- lapply(1:10, function(x){
lmer_fun2(X = x, DV = dv, moder = Moder, form = Form, Groups = groups, PSW.Weights = psw.weights)})
pred.fx <- fx_pred_fun(mods, Form)
save(pred.fx, file = sprintf("%s/results/%s/pred/%s_%s_%s_%s.RData",
data_path, Form, psw.weights, dv, groups, Moder))
# p <- pool_fun(mods)
# save(p, file = sprintf("%s/results/%s/pooled/%s_%s_%s_%s_%s.RData",
# data_path, Form, psw.weights, dv, groups, Moder))
# p_re <- pooled_re_fun(mods)
# save(p_re, file = sprintf("%s/results/%s/pooled_re/%s_%s_%s_%s.RData",
# data_path, Form, psw.weights, dv, groups, Moder))
rm(list = c("mods", "p", "p_re", "pred"))
gc()
}
# short function for running models
lmer_fun2 <- function(X, DV, moder, form, Groups, PSW.Weights){
print(paste(DV, moder, form, sep = " "))
load(sprintf("%s/data/model_data/%s/%s_chain.%s_%s.RData", data_path, form, DV, X, Groups))
if(PSW.Weights == "weighted"){pw <- d$psw.weights}
else{pw <- rep(1, nrow(d))}
md <- moder_wide %>%
select(PROC_CID, one_of(moder)) %>%
setNames(c("PROC_CID", "modvalue"))
if(class(md$modvalue) != "factor"){md <- md %>% mutate(modvalue = as.numeric(scale(modvalue, scale = F)))}
mage <- d %>% select(PROC_CID, cont_age, cont_year) %>% distinct() %>% summarize_at(vars(-PROC_CID), mean, na.rm = T)
d <- d %>%
left_join(md) %>%
group_by(PROC_CID) %>%
mutate(age_c = age_baseline - mage$cont_age,
wave0 = year - mage$cont_year) %>%
ungroup() %>%
mutate(groups = fct_relevel(groups, "None"))
cl <- lmerControl(check.conv.grad = .makeCC("ignore", tol = 2e-2, relTol = NULL),
check.conv.singular = .makeCC(action = "ignore", tol = 1e-2),
check.conv.hess = .makeCC(action = "ignore", tol = 1e-3),
check.scaleX = c("ignore"))
y <- if(form == "linear"){
lmer(value ~ T1_value + age0*groups*modvalue + (age0 | PROC_CID)
, data = d, weights = pw, control = cl)
} else{
lmer(value ~ wave0 + age_c + tv_group + groups + modvalue + wave0:groups +
age_c:groups + wave0:tv_group + wave0:modvalue + tv_group:modvalue +
modvalue:groups + modvalue:wave0:groups + (wave0 | PROC_CID)
, data = d, weights = pw, control = cl)
}
save(y, file = sprintf("%s/results/%s/models/%s_%s_chain.%s_%s_%s.RData", data_path, form, PSW.Weights, DV, X, Groups, moder))
rm(d)
gc()
return(y)
}7.1.2 Prediction Function
fx_pred_fun <- function(fits, m){
if(m == "linear"){
x <- lapply(fits, function(x) {
cl <- class(x@frame$modvalue)
if(cl %in% c("numeric", "integer", "matrix")){
msd <- x@frame %>%
select(modvalue) %>%
as_tibble() %>%
summarize_all(lst(mean, sd))
pred.fx <- crossing(
age0 = seq(-8,20,.2)/5
, groups = c("None", "CommServ", "Jail")
, modvalue = with(msd, c(mean-sd, mean,mean+sd))
, T1_value = mean(x@frame$T1_value)
)
} else {
pred.fx <- crossing(
age0 = seq(-8,20,.2)/5
, groups = c("None", "CommServ", "Jail")
, modvalue = factor(levels(x@frame$modvalue))
, T1_value = mean(x@frame$T1_value)
)
}
pred.fx <- pred.fx %>%
mutate(pred = predict(x
, newdata = pred.fx
, re.form = NA))
if(cl %in% c("numeric", "integer", "matrix")){
pred.fx <- pred.fx %>% mutate(modvalue = mapvalues(modvalue, unique(sort(modvalue)), c("-1 SD", "M", "+1 SD")))
}
return(pred.fx)
})
res <- tibble(pred = x) %>%
unnest(pred) %>%
group_by(age0, groups, modvalue, T1_value) %>%
summarise_all(mean, na.rm = T) %>%
ungroup()
} else {
x <- lapply(fits, function(x) {
cl <- class(x@frame$modvalue)
if(cl %in% c("numeric", "integer", "matrix")){
msd <- x@frame %>%
select(modvalue) %>%
as_tibble() %>%
summarize_all(lst(mean, sd))
pred.fx <- crossing(
wave0 = seq(-10, 10, .001)
, groups = factor(levels(x@frame$groups))
, modvalue = with(msd, c(mean-sd, mean,mean+sd))
, age_c = 0
, tv_group = c(0,1)
) %>%
filter((groups == "None" & tv_group == 0) |
(groups != "None" & tv_group == 0 & wave0 < 0) |
(groups != "None" & tv_group == 1 & wave0 > 0))
} else {
pred.fx <- crossing(
wave0 = seq(-10, 10, .001)
, groups = factor(levels(x@frame$groups))
, modvalue = factor(levels(x@frame$modvalue))
, age_c = 0
, tv_group = c(0,1)
) %>%
filter((groups == "None" & tv_group == 0) |
(groups != "None" & tv_group == 0 & wave0 < 0) |
(groups != "None" & tv_group == 1 & wave0 > 0))
}
pred.fx <- pred.fx %>%
mutate(pred = predict(x
, newdata = pred.fx
, re.form = NA)) %>%
mutate(pred = ifelse(wave0 == 0, NA, pred))
if(cl %in% c("numeric", "integer", "matrix")){
pred.fx <- pred.fx %>% mutate(modvalue = mapvalues(modvalue, unique(sort(modvalue)), c("-1 SD", "M", "+1 SD")))
}
return(pred.fx)
})
res <- tibble(pred = x) %>%
unnest(pred) %>%
group_by(wave0, groups, modvalue, age_c, tv_group) %>%
summarise_all(mean, na.rm = T) %>%
ungroup()
}
return(res)
}
loadRData <- function(path, obj){
#loads an RData file, and returns it
load(path)
get(ls()[grepl(obj, ls())])
}7.1.3 Hypothesis Function
moder_hypoth_fun <- function(fits, m){
y <- lapply(fits, function(x){
cl <- class(x@frame$modvalue)
if(m == "linear"){
if(cl == "factor"){
K <- c("age0 = 0", # moder level 0, no contact
"age0 + age0:groupsJail = 0", # moder level 0, jail
"age0 + age0:groupsCommServ = 0", # moder level 0, commserv
"age0 + age0:modvalue1 = 0", # moder level 1, no contact
"age0 + age0:modvalue1 + age0:groupsJail:modvalue1 = 0", # moder level 1, jail
"age0 + age0:modvalue1 + age0:groupsCommServ:modvalue1 = 0") # moder level 1, jail
grnm <- c("None_modvalue0", "Jail_modvalue0", "CommServ_modvalue0", "None_modvalue1",
"Jail_modvalue1", "CommServ_modvalue1")
if(length(levels(x@frame$modvalue)) == 3){
K <- c(K,
"age0 + age0:modvalue2 = 0", # moder level 1, no contact
"age0 + age0:modvalue2 + age0:groupsJail:modvalue2 = 0", # moder level 1, jail
"age0 + age0:modvalue2 + age0:groupsCommServ:modvalue2 = 0") # moder level 1, jail
grnm <- c(grnm, "None_modvalue2", "Jail_modvalue2", "CommServ_modvalue2")
}
} else {
msd <- x@frame %>%
select(modvalue) %>%
as_tibble() %>%
summarize_all(lst(mean, sd)) %>%
mutate(lower = mean - sd, upper = mean + sd)
K <- rbind("None_modvalue_M" = c(0,0,1,0,0,0,0,0,msd$mean,0,0,0,0),
"Jail_modvalue_M" = c(0,0,1,0,0,1,0,0,msd$mean,0,0,msd$mean,0),
"CommServ_modvalue_M" = c(0,0,1,0,0,0,1,0,msd$mean,0,0,0,msd$mean),
"None_modvalue_lower" = c(0,0,1,0,0,0,0,0,msd$lower,0,0,0,0),
"Jail_modvalue_lower" = c(0,0,1,0,0,1,0,0,msd$lower,0,0,msd$lower,0),
"CommServ_modvalue_lower" = c(0,0,1,0,0,0,1,0,msd$lower,0,0,0,msd$lower),
"None_modvalue_upper" = c(0,0,1,0,0,0,0,0,msd$upper,0,0,0,0),
"Jail_modvalue_upper" = c(0,0,1,0,0,1,0,0,msd$upper,0,0,msd$upper,0),
"CommServ_modvalue_upper" = c(0,0,1,0,0,0,1,0,msd$upper,0,0,0,msd$upper)
)
grnm <- rownames(K)
}
} else {
gr <- names(fixef(x))[grepl('groups', names(fixef(x)))][1]
if(cl == "factor"){
K <- if(length(levels(x@frame$modvalue)) == 2){
rbind("Pre_modvalue0" = c(1,0,0, 0,1,0, 0,0,0, 0,0,0, 0),
"Post_modvalue0" = c(1,0,0, 1,1,0, 0,0,0, 0,0,0, 0),
"Pre_modvalue1" = c(1,0,0, 0,1,1, 0,0,0, 0,0,1, 0),
"Post_modvalue1" = c(1,0,0, 1,1,1, 0,0,0, 0,1,1, 0))
} else {
rbind("Pre_modvalue0" = c(1,0,0, 0,1,0, 0,0,0, 0,0,0, 0,0,0, 0,0,0),
"Post_modvalue0" = c(1,0,0, 1,1,0, 0,0,0, 0,0,0, 0,0,0, 0,0,0),
"Pre_modvalue1" = c(1,0,0, 0,1,1, 0,0,0, 0,0,0, 0,0,1, 0,0,0),
"Post_modvalue1" = c(1,0,0, 1,1,1, 0,0,0, 0,0,0, 1,0,1, 0,0,0),
"Pre_modvalue2" = c(1,0,0, 0,1,0, 1,0,0, 0,0,0, 0,0,0, 1,0,0),
"Post_modvalue2" = c(1,0,0, 1,1,0, 1,0,0, 0,0,0, 0,1,0, 1,0,0))
}
grnm <- rownames(K)
} else {
msd <- x@frame %>%
select(modvalue) %>%
as_tibble() %>%
summarize_all(lst(mean, sd)) %>%
mutate(lower = mean - sd, upper = mean + sd)
K <- rbind("Pre_modvalue_M" = c(1,0,0, 0,1,msd$mean, 0,0,0, 0,0,msd$mean, 0),
"Post_modvalue_M" = c(1,0,0, 1,1,msd$mean, 0,0,0, 0,msd$mean,msd$mean, 0),
"Pre_modvalue_lower" = c(1,0,0, 0,1,msd$lower, 0,0,0, 0,0,msd$lower, 0),
"Post_modvalue_lower" = c(1,0,0, 1,1,msd$lower, 0,0,0, 0,msd$lower,msd$lower, 0),
"Pre_modvalue_upper" = c(1,0,0, 0,1,msd$upper, 0,0,0, 0,0,msd$upper, 0),
"Post_modvalue_upper" = c(1,0,0, 1,1,msd$upper, 0,0,0, 0,msd$upper,msd$upper, 0)
)
grnm <- rownames(K)
}
}
h <- confint(multcomp::glht(x, linfct = K))$confint %>%
data.frame() %>%
mutate(term = grnm)
})
res <- tibble(hypoth = y) %>%
unnest(hypoth) %>%
group_by(term) %>%
summarise_all(mean, na.rm = T) %>%
ungroup()
return(res)
}7.2 Run and Pool Models
nested.moder.mods <- crossing(
pred = unique(traits$old),
mod = c("linear", "spline"),
moder = unique(moder_long$moditem),
groups = c("Jail", "CommServ", "threegroup"),
weighted = c("weighted", "unweighted")
) %>%
filter(!(mod == "linear" & groups %in% c("Jail", "CommServ")) &
!(mod == "spline" & groups == "threegroup")) %>%
filter(!(moder %in% small_size)) %>%
filter(moder %in% c("SclSat", "SclChal", "EduRemed")) %>%
# filter(weighted == "weighted" &
# ((pred %in% c("SensSeek", "Impls") & groups == "CommServ" & moder == "DemRace") |
# (pred == "CESD" & groups == "CommServ" & moder %in% c("DemRace", "SclChal", "SclSat", "ParDvrcd", "EduRemed", "DSSpank", "DemCSex")))) %>%
# select(pred, moder, mod, groups, weighted ) %>%
# write.table(., file = sprintf("%s/scripts/cluster/moder_args.txt", data_path), col.names = F, row.names = F)
mutate(models = pmap(list(pred, moder, mod, groups, weighted) , mod_fun2))
gc()7.3 Load Results
# clear up ram
loadRData <- function(fileName, form, type, obj){
#loads an RData file, and returns it
path <- sprintf("%s/results/%s/%s/%s", local_path, form, type, fileName)
load(path)
get(ls()[grepl(obj, ls())])
}
nested.moder.mods <- tibble(mod = c("spline", "linear")) %>%
mutate(file = map(mod, ~list.files(sprintf("%s/results/%s/pooled", local_path,.)))) %>%
unnest(file) %>%
separate(file, c("weighted", "pred", "groups", "moder"), sep = "_", remove = F) %>%
mutate(moder = str_remove_all(moder, ".RData")) %>%
filter(moder != "none") %>%
mutate(pool = map2(file, mod, ~loadRData(.x, .y, "pooled", "p")),
pool_re = map2(file, mod, ~loadRData(.x, .y, "pooled_re", "p_re")),
pred.fx = map2(file, mod, ~loadRData(.x, .y, "pred", "pred.fx")),
hypoth = map2(file, mod, ~loadRData(.x, .y, "hypoth", "hypoth")))7.4 Table Results
mod.tab <- nested.moder.mods %>%
select(mod, weighted:pool) %>%
unnest(pool) %>%
# get rid of the random effects for now
filter(type == "fixeff") # %>%
# mutate(CI = sprintf("[%.2f, %.2f]", Estimate - CI, Estimate + CI))
results.tab <- mod.tab %>%
#### fixed effects ####
mutate(type = ifelse(grepl("modvalue", term) == T, "Moderator", type)) %>%
# format the confidence intervals and estimate for LaTeX
mutate(CI = sprintf("[%.2f, %.2f]", CI.lower, CI.upper),
CI = ifelse(p < .05, sprintf("<strong>%s</strong>", CI), CI),
Estimate = round(Estimate, 2),
Estimate = ifelse(p < .05, sprintf("<strong>%s</strong>", Estimate), Estimate)) %>%
# get rid of t, p, and the empty group variable
select(-t, -p, -CI.lower, -CI.upper) %>%
rename(b = Estimate) %>%
# change b & CI to long so we can make the outcomes, rather than coefficients, wide
gather(key = est, value = value, b, CI)7.4.1 Linear Models
moders2 <- tribble(
~new_Moder, ~Moder, ~old, ~new,
"Gender", "DemCSex", "0", "Male",
"Gender", "DemCSex", "1", "Female",
"Race", "DemRace", "0", "White",
"Race", "DemRace", "1", "Black",
"Race", "DemRace", "2", "Other",
"Corporal Punishment", "DSSpank", "0", "No",
"Corporal Punishment", "DSSpank", "1", "Yes",
"Remedial Education", "EduRemed", "-1 SD", "-1 SD",
"Remedial Education", "EduRemed", "M", "M",
"Remedial Education", "EduRemed", "+1 SD", "+1 SD",
"Parent Died", "HHParDead", "0", "None",
"Parent Died", "HHParDead", "1", "One",
"Parent Died", "HHParDead", "2", "Two",
"Welfare", "MWelfare", "0", "No",
"Welfare", "MWelfare", "1", "Yes",
"Parental Divorce", "ParDvrcd", "0", "No",
"Parental Divorce", "ParDvrcd", "1", "Yes",
"School Challenging", "SclChal", "-1 SD", "-1 SD",
"School Challenging", "SclChal", "M", "M",
"School Challenging", "SclChal", "+1 SD", "+1 SD",
"School Satisfaction", "SclSat", "-1 SD", "-1 SD",
"School Satisfaction", "SclSat", "M", "M",
"School Satisfaction", "SclSat", "+1 SD", "+1 SD"
) %>%
mutate(comb = paste0(Moder, old))
moders3 <- tribble(
~new_Moder, ~old, ~new,
"Gender", "DemCSex1", "Male v Female",
"Race", "DemRace1", "White v Black",
"Race", "DemRace2", "White v Other",
"Corporal Punishment", "DSSpank1", "No v Yes",
"Parent Died", "HHParDead1", "None v One",
"Parent Died", "HHParDead2", "None v Two",
"Welfare", "MWelfare1", "No v Yes",
"Parental Divorce", "ParDvrcd1", "No v Yes",
"School Satisfaction", "SclSat", "Continuous",
"School Challenging", "SclChal", "Continuous",
"Remedial Education", "EduRemed", "Years"
)
lin_moder_tab_fun <- function(d, weight){
options(papaja.na_string = "")
sup_num <- if(weight == "weighted") 4 else paste0("S",sup_tables)
cap <- sprintf("<strong>Table %s</strong><br><em>Moderators of Linear Slopes between Groups for %s Models", sup_num, weight)
tab <- d %>%
mutate(groups = factor(groups, levels = c("Jail", "CommServ")),
moder = factor(moder, levels = unique(moders2$new_Moder)),
term = factor(term, levels = unique(d$term))) %>%
arrange(groups, moder, term) %>%
select(-groups) %>%
kable(., "html", booktabs = T, escape = F,
col.names = c("", "Group", rep("b [CI]", times = 4)),
caption = cap,
align = c("l", "l", rep("c", 4))) %>%
kable_styling(full_width = F) %>%
add_header_above(c(" " = 2, "Sensation Seeking" = 1, "Impulsivity" = 1, "Depressive Symptoms" = 1, "Self Esteem" = 1)) %>%
kableExtra::group_rows("None v Incarceration", 1, 11) %>%
kableExtra::group_rows("None v Community Service", 12, 22) %>%
collapse_rows(1, valign = "top")
save_kable(tab, file = sprintf("%s/results/tables/Table_%s_%s_bw_person_moder_slopes.html", local_path, sup_num, weight))
if(weight == "unweighted") sup_tables <<- sup_tables + 1
return(tab)
}
lin_moder_tab <- results.tab %>%
filter((grepl("modvalue", term) & grepl("^age0:gr", term)) & mod == "linear") %>%
separate(term, into = c("t1", "t2", "t3"), sep = ":") %>%
mutate(t3 = str_replace_all(t3, "modvalue", moder),
t2 = str_remove_all(t2, "groups")) %>%
arrange(pred, t2, est) %>%
select(-d) %>%
pivot_wider(names_from = "est", values_from = "value") %>%
mutate(moder = mapvalues(moder, moders2$Moder, moders2$new_Moder),
t3 = mapvalues(t3, moders3$old, moders3$new),
estimate = sprintf("%s<br>%s", b, CI)) %>%
select(-b, -CI) %>%
pivot_wider(names_from = "pred", values_from = "estimate") %>%
select(groups = t2, weighted, moder, term = t3, contains(traits$old))
lin_moder_kable <- lin_moder_tab %>%
group_by(weighted) %>%
nest() %>%
ungroup() %>%
mutate(tab = map2(data, weighted, lin_moder_tab_fun))7.4.1.1 Weighted
| Group | b [CI] | b [CI] | b [CI] | b [CI] | |
|---|---|---|---|---|---|
| None v Incarceration | |||||
| Gender | Male v Female |
-0.04 [-0.10, 0.02] |
0.05 [-0.01, 0.10] |
-0.09 [-0.15, -0.03] |
0.02 [-0.03, 0.07] |
| Race | White v Black |
-0.01 [-0.07, 0.06] |
0.02 [-0.04, 0.08] |
0.1 [0.04, 0.16] |
0.01 [-0.04, 0.06] |
| White v Other |
-0.12 [-0.18, -0.05] |
-0.03 [-0.09, 0.04] |
0.02 [-0.04, 0.08] |
0.07 [0.02, 0.13] |
|
| Corporal Punishment | No v Yes |
-0.01 [-0.07, 0.05] |
0.03 [-0.03, 0.09] |
0 [-0.06, 0.05] |
-0.01 [-0.06, 0.04] |
| Remedial Education | Years |
-0.01 [-0.07, 0.04] |
0.03 [-0.03, 0.08] |
0 [-0.06, 0.05] |
0.01 [-0.04, 0.05] |
| Parent Died | None v One |
0.06 [0.00, 0.12] |
-0.03 [-0.08, 0.03] |
0.04 [-0.01, 0.09] |
-0.01 [-0.05, 0.03] |
| None v Two |
0.04 [-0.09, 0.18] |
-0.04 [-0.17, 0.09] |
-0.01 [-0.14, 0.13] |
0.04 [-0.07, 0.15] |
|
| Welfare | No v Yes |
-0.07 [-0.13, -0.01] |
-0.13 [-0.19, -0.08] |
0.07 [0.02, 0.13] |
0.09 [0.04, 0.14] |
| Parental Divorce | No v Yes |
0.02 [-0.04, 0.08] |
0.02 [-0.03, 0.07] |
0.01 [-0.05, 0.06] |
-0.01 [-0.05, 0.03] |
| School Challenging | Continuous |
0.01 [-0.07, 0.09] |
0.01 [-0.06, 0.09] |
0.02 [-0.06, 0.09] |
-0.06 [-0.11, -0.00] |
| School Satisfaction | Continuous |
-0.01 [-0.07, 0.05] |
0.02 [-0.04, 0.08] |
-0.01 [-0.07, 0.05] |
0 [-0.04, 0.05] |
| None v Community Service | |||||
| Gender | Male v Female |
0.01 [-0.06, 0.08] |
0.01 [-0.06, 0.07] |
-0.02 [-0.08, 0.04] |
0 [-0.06, 0.06] |
| Race | White v Black |
-0.01 [-0.08, 0.07] |
0.07 [-0.00, 0.15] |
0 [-0.08, 0.07] |
0.02 [-0.04, 0.09] |
| White v Other |
-0.1 [-0.19, -0.02] |
-0.07 [-0.15, 0.01] |
0.03 [-0.05, 0.10] |
0.02 [-0.05, 0.09] |
|
| Corporal Punishment | No v Yes |
-0.03 [-0.11, 0.05] |
0 [-0.08, 0.08] |
0.09 [0.01, 0.17] |
0.01 [-0.05, 0.07] |
| Remedial Education | Years |
-0.05 [-0.12, 0.02] |
0.01 [-0.06, 0.07] |
-0.05 [-0.12, 0.02] |
0.06 [0.00, 0.12] |
| Parent Died | None v One |
-0.01 [-0.08, 0.06] |
-0.02 [-0.08, 0.05] |
-0.07 [-0.13, -0.00] |
0.09 [0.03, 0.15] |
| None v Two |
0.2 [0.03, 0.36] |
0.15 [0.00, 0.31] |
0 [-0.15, 0.15] |
-0.03 [-0.16, 0.11] |
|
| Welfare | No v Yes |
0 [-0.07, 0.07] |
-0.04 [-0.11, 0.03] |
0.03 [-0.04, 0.09] |
0.03 [-0.04, 0.09] |
| Parental Divorce | No v Yes |
-0.05 [-0.12, 0.02] |
-0.01 [-0.08, 0.06] |
0.04 [-0.03, 0.10] |
0.03 [-0.02, 0.09] |
| School Challenging | Continuous |
-0.04 [-0.13, 0.04] |
-0.11 [-0.20, -0.02] |
0.03 [-0.04, 0.11] |
-0.08 [-0.16, -0.01] |
| School Satisfaction | Continuous |
-0.08 [-0.16, 0.01] |
-0.01 [-0.09, 0.06] |
0.1 [0.03, 0.18] |
-0.03 [-0.10, 0.04] |
7.4.1.2 Unweighted
| Group | b [CI] | b [CI] | b [CI] | b [CI] | |
|---|---|---|---|---|---|
| None v Incarceration | |||||
| Gender | Male v Female |
-0.05 [-0.13, 0.03] |
0.04 [-0.04, 0.11] |
-0.03 [-0.10, 0.03] |
0.03 [-0.03, 0.09] |
| Race | White v Black |
-0.01 [-0.09, 0.06] |
0.02 [-0.05, 0.09] |
0.05 [-0.01, 0.11] |
0.01 [-0.04, 0.07] |
| White v Other |
-0.11 [-0.19, -0.03] |
-0.04 [-0.11, 0.04] |
0.01 [-0.06, 0.07] |
0.04 [-0.02, 0.10] |
|
| Corporal Punishment | No v Yes |
-0.02 [-0.09, 0.05] |
0.01 [-0.06, 0.07] |
-0.02 [-0.08, 0.04] |
-0.01 [-0.06, 0.04] |
| Remedial Education | Years |
-0.01 [-0.08, 0.05] |
0.02 [-0.04, 0.08] |
-0.01 [-0.07, 0.04] |
-0.03 [-0.08, 0.02] |
| Parent Died | None v One |
0.07 [0.01, 0.14] |
-0.04 [-0.10, 0.02] |
-0.02 [-0.08, 0.03] |
-0.01 [-0.05, 0.04] |
| None v Two |
0.06 [-0.08, 0.20] |
-0.04 [-0.17, 0.10] |
-0.06 [-0.18, 0.07] |
0.05 [-0.05, 0.15] |
|
| Welfare | No v Yes |
-0.03 [-0.11, 0.05] |
-0.1 [-0.17, -0.03] |
0.04 [-0.02, 0.11] |
0.05 [-0.01, 0.11] |
| Parental Divorce | No v Yes |
0.03 [-0.04, 0.09] |
0.04 [-0.02, 0.10] |
-0.01 [-0.06, 0.05] |
-0.01 [-0.06, 0.04] |
| School Challenging | Continuous |
0.01 [-0.07, 0.10] |
0.05 [-0.03, 0.13] |
0 [-0.07, 0.08] |
-0.05 [-0.11, 0.00] |
| School Satisfaction | Continuous |
0 [-0.06, 0.07] |
0.01 [-0.05, 0.07] |
-0.02 [-0.07, 0.04] |
0.01 [-0.04, 0.06] |
| None v Community Service | |||||
| Gender | Male v Female |
0.01 [-0.08, 0.10] |
0.02 [-0.07, 0.10] |
-0.02 [-0.10, 0.05] |
-0.01 [-0.07, 0.06] |
| Race | White v Black |
0.09 [-0.01, 0.20] |
0.06 [-0.04, 0.15] |
-0.03 [-0.11, 0.06] |
0.01 [-0.06, 0.09] |
| White v Other |
-0.06 [-0.17, 0.04] |
-0.06 [-0.16, 0.04] |
-0.01 [-0.10, 0.08] |
0.02 [-0.06, 0.10] |
|
| Corporal Punishment | No v Yes |
0.01 [-0.09, 0.11] |
0 [-0.10, 0.09] |
0.02 [-0.06, 0.11] |
0.03 [-0.04, 0.09] |
| Remedial Education | Years |
-0.02 [-0.11, 0.07] |
0.01 [-0.07, 0.10] |
-0.05 [-0.12, 0.03] |
0.05 [-0.02, 0.12] |
| Parent Died | None v One |
0.04 [-0.05, 0.13] |
0 [-0.09, 0.08] |
-0.09 [-0.17, -0.01] |
0.07 [-0.01, 0.14] |
| None v Two |
0.21 [0.01, 0.41] |
0.08 [-0.10, 0.27] |
-0.09 [-0.26, 0.08] |
0.04 [-0.11, 0.19] |
|
| Welfare | No v Yes |
0.05 [-0.05, 0.15] |
0 [-0.10, 0.09] |
-0.02 [-0.10, 0.07] |
0.02 [-0.06, 0.10] |
| Parental Divorce | No v Yes |
-0.03 [-0.12, 0.06] |
-0.02 [-0.11, 0.06] |
0.04 [-0.04, 0.12] |
0.02 [-0.05, 0.09] |
| School Challenging | Continuous |
-0.06 [-0.17, 0.05] |
-0.07 [-0.17, 0.03] |
0.02 [-0.07, 0.11] |
-0.06 [-0.14, 0.02] |
| School Satisfaction | Continuous |
-0.03 [-0.13, 0.07] |
-0.01 [-0.11, 0.09] |
0.11 [0.02, 0.19] |
-0.02 [-0.10, 0.05] |
7.4.1.3 Simple Effects
lin_moder_sf_tab <- function(d, weighted){
spanh <- c(2, rep(1,12))
names(spanh) <- c(" ", rep(c("Neither", "Incarceration", "Community\nService"), times = 4))
tab <- d %>%
kable(.
, "html"
, escape = F
, col.names = c("Moderator", "", rep("b [CI]", 12))
, align = c("r", "r", rep("c", 12))) %>%
kable_styling(full_width = F) %>%
add_header_above(spanh) %>%
add_header_above(c(" " = 2, "Sensation Seeking" = 3, "Impulsivity" = 3,
"Depressive Symptoms" = 3, "Self-Esteem" = 3)) %>%
collapse_rows(1, valign = "top")
save_kable(tab, file = sprintf("%s/results/tables/Table_S%s_%s_bw_person_moder_simple_effects.html", local_path, sup_tables, weighted))
sup_tables <<- sup_tables + 1
return(tab)
}
levs <- paste(rep(traits$old, each = 3), c("None", "Jail", "CommServ"), sep = "_")
lin_moder_sf <- nested.moder.mods %>%
filter(mod == "linear") %>%
select(mod, weighted:moder, hypoth) %>%
unnest(hypoth) %>%
separate(term, c("term2", "term", "modval"), sep = "_") %>%
mutate(modval = ifelse(is.na(modval), str_remove_all(term, "[A-Z a-z]"), modval),
modval = mapvalues(modval, c("lower", "upper"), c("-1 SD", "+1 SD")),
modval = paste0(moder, modval),
modval = mapvalues(modval, moders2$comb, moders2$new)) %>%
mutate(Estimate = sprintf("%.2f<br>[%.2f, %.2f]", Estimate, lwr, upr)) %>%
select(weighted:Estimate) %>%
left_join(
nested.moder.mods %>%
filter(mod == "linear") %>%
select(mod, weighted:moder, pool) %>%
unnest(pool) %>%
filter(grepl("age0:groupsJail:", term)|grepl("age0:groupsCommServ:", term)) %>%
mutate(term = str_remove_all(term, "age0:groups"),
sig = ifelse(sign(CI.lower)==sign(CI.upper), "sig", "ns")) %>%
separate(term, c("term2", "term"), sep = ":") %>%
select(weighted:pred, moder, term, term2, sig)
) %>%
mutate(Estimate = ifelse(sig == "sig" & !is.na(sig),
sprintf("<strong>%s</strong>", Estimate), Estimate),
moder = factor(moder, moders2$Moder, moders2$new_Moder)) %>%
select(-sig, -term) %>%
pivot_wider(names_from = c("pred", "term2")
, values_from = "Estimate"
, names_sep = "_") %>%
select(weighted:modval, one_of(levs), -groups) %>%
group_by(weighted) %>%
nest() %>%
ungroup() %>%
mutate(tab = map2(data, weighted, lin_moder_sf_tab))| Moderator | b [CI] | b [CI] | b [CI] | b [CI] | b [CI] | b [CI] | b [CI] | b [CI] | b [CI] | b [CI] | b [CI] | b [CI] | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Gender | Male |
0.03 [0.02, 0.05] |
0.04 [0.01, 0.08] |
0.06 [0.01, 0.12] |
-0.10 [-0.11, -0.08] |
-0.07 [-0.11, -0.04] |
-0.08 [-0.13, -0.03] |
-0.01 [-0.03, 0.00] |
-0.03 [-0.07, 0.00] |
0.01 [-0.04, 0.06] |
0.04 [0.03, 0.06] |
0.05 [0.02, 0.07] |
0.07 [0.03, 0.11] |
| Female |
0.01 [-0.01, 0.03] |
-0.03 [-0.10, 0.05] |
0.03 [-0.06, 0.11] |
-0.11 [-0.13, -0.10] |
-0.07 [-0.14, 0.01] |
-0.11 [-0.19, -0.02] |
-0.02 [-0.04, -0.00] |
-0.11 [-0.18, -0.04] |
-0.04 [-0.12, 0.04] |
0.05 [0.04, 0.06] |
0.07 [0.01, 0.13] |
0.05 [-0.02, 0.12] |
|
| Race | White |
0.02 [-0.00, 0.04] |
0.05 [-0.00, 0.11] |
0.07 [0.01, 0.14] |
-0.09 [-0.11, -0.07] |
-0.05 [-0.11, -0.00] |
-0.07 [-0.13, -0.01] |
-0.02 [-0.04, -0.00] |
-0.10 [-0.15, -0.05] |
-0.01 [-0.07, 0.05] |
0.04 [0.02, 0.05] |
0.02 [-0.02, 0.07] |
0.05 [-0.00, 0.10] |
| Black |
0.02 [-0.00, 0.04] |
0.01 [-0.07, 0.09] |
0.01 [-0.09, 0.12] |
-0.12 [-0.14, -0.10] |
-0.10 [-0.18, -0.02] |
-0.05 [-0.15, 0.05] |
0.00 [-0.02, 0.02] |
0.10 [0.02, 0.17] |
-0.00 [-0.10, 0.09] |
0.04 [0.03, 0.06] |
0.05 [-0.01, 0.12] |
0.06 [-0.02, 0.15] |
|
| Other |
0.06 [0.04, 0.09] |
-0.05 [-0.14, 0.03] |
-0.04 [-0.15, 0.07] |
-0.11 [-0.13, -0.08] |
-0.14 [-0.22, -0.05] |
-0.18 [-0.28, -0.08] |
-0.04 [-0.06, -0.01] |
-0.02 [-0.10, 0.06] |
-0.01 [-0.10, 0.08] |
0.07 [0.04, 0.09] |
0.14 [0.07, 0.20] |
0.09 [0.00, 0.17] |
|
| Corporal Punishment | No |
-0.01 [-0.03, 0.01] |
-0.01 [-0.07, 0.06] |
0.04 [-0.05, 0.13] |
-0.10 [-0.12, -0.08] |
-0.08 [-0.14, -0.02] |
-0.07 [-0.15, 0.01] |
-0.01 [-0.03, 0.01] |
-0.06 [-0.11, 0.00] |
-0.07 [-0.14, 0.01] |
0.04 [0.02, 0.05] |
0.05 [0.00, 0.10] |
0.05 [-0.01, 0.12] |
| Yes |
0.04 [0.03, 0.06] |
0.04 [-0.04, 0.11] |
0.01 [-0.09, 0.11] |
-0.11 [-0.12, -0.09] |
-0.08 [-0.15, -0.00] |
-0.11 [-0.21, -0.02] |
-0.02 [-0.03, -0.00] |
-0.02 [-0.09, 0.05] |
0.07 [-0.02, 0.16] |
0.05 [0.04, 0.06] |
0.04 [-0.02, 0.10] |
0.06 [-0.01, 0.14] |
|
| Remedial Education | -1 SD |
0.05 [0.03, 0.07] |
0.03 [-0.03, 0.08] |
0.06 [-0.01, 0.14] |
-0.12 [-0.14, -0.10] |
-0.08 [-0.13, -0.03] |
-0.10 [-0.17, -0.02] |
-0.04 [-0.05, -0.02] |
-0.02 [-0.08, 0.03] |
-0.04 [-0.10, 0.02] |
0.05 [0.03, 0.06] |
-0.04 [-0.08, 0.01] |
0.02 [-0.03, 0.08] |
| M |
0.06 [0.05, 0.07] |
0.05 [0.01, 0.08] |
0.04 [-0.01, 0.08] |
-0.11 [-0.12, -0.10] |
-0.06 [-0.09, -0.03] |
-0.09 [-0.13, -0.04] |
-0.03 [-0.04, -0.02] |
-0.01 [-0.04, 0.03] |
-0.06 [-0.10, -0.02] |
0.05 [0.04, 0.06] |
-0.02 [-0.05, 0.01] |
0.05 [0.01, 0.08] |
|
| +1 SD |
0.07 [0.05, 0.09] |
0.06 [0.01, 0.12] |
0.01 [-0.07, 0.08] |
-0.11 [-0.13, -0.09] |
-0.03 [-0.09, 0.02] |
-0.08 [-0.15, -0.01] |
-0.02 [-0.03, -0.00] |
0.01 [-0.04, 0.06] |
-0.09 [-0.15, -0.03] |
0.06 [0.04, 0.07] |
-0.00 [-0.05, 0.04] |
0.07 [0.02, 0.12] |
|
| Parent Died | None |
0.05 [0.03, 0.07] |
0.02 [-0.04, 0.07] |
0.07 [0.01, 0.14] |
-0.10 [-0.12, -0.08] |
-0.05 [-0.10, 0.00] |
-0.07 [-0.14, -0.01] |
-0.02 [-0.03, 0.00] |
-0.08 [-0.13, -0.03] |
0.03 [-0.02, 0.09] |
0.05 [0.03, 0.06] |
0.06 [0.02, 0.10] |
0.02 [-0.04, 0.08] |
| One |
0.02 [0.00, 0.04] |
0.08 [0.00, 0.15] |
0.01 [-0.08, 0.10] |
-0.11 [-0.13, -0.09] |
-0.13 [-0.20, -0.06] |
-0.13 [-0.22, -0.04] |
-0.03 [-0.04, -0.01] |
0.01 [-0.05, 0.08] |
-0.09 [-0.17, -0.01] |
0.05 [0.04, 0.06] |
0.04 [-0.02, 0.10] |
0.14 [0.06, 0.21] |
|
| Two |
0.02 [-0.05, 0.08] |
0.06 [-0.11, 0.24] |
0.22 [0.00, 0.43] |
-0.12 [-0.18, -0.05] |
-0.16 [-0.32, 0.01] |
0.04 [-0.17, 0.24] |
-0.02 [-0.08, 0.05] |
-0.02 [-0.19, 0.14] |
-0.01 [-0.20, 0.17] |
0.03 [-0.02, 0.07] |
0.07 [-0.06, 0.20] |
0.00 [-0.16, 0.16] |
|
| Welfare | No |
0.05 [0.03, 0.07] |
0.11 [0.05, 0.17] |
0.08 [0.01, 0.16] |
-0.12 [-0.13, -0.10] |
0.02 [-0.04, 0.07] |
-0.06 [-0.13, 0.01] |
-0.03 [-0.05, -0.01] |
-0.12 [-0.17, -0.07] |
-0.03 [-0.09, 0.03] |
0.05 [0.03, 0.07] |
-0.01 [-0.06, 0.04] |
0.05 [-0.01, 0.11] |
| Yes |
0.02 [0.00, 0.03] |
-0.06 [-0.13, 0.02] |
0.02 [-0.07, 0.11] |
-0.11 [-0.12, -0.09] |
-0.24 [-0.31, -0.17] |
-0.15 [-0.24, -0.06] |
-0.01 [-0.03, -0.00] |
0.06 [-0.01, 0.12] |
0.01 [-0.07, 0.09] |
0.05 [0.04, 0.06] |
0.13 [0.08, 0.19] |
0.08 [0.00, 0.15] |
|
| Parental Divorce | No |
0.03 [0.01, 0.04] |
0.03 [-0.01, 0.07] |
0.08 [0.03, 0.13] |
-0.10 [-0.12, -0.09] |
-0.07 [-0.11, -0.03] |
-0.07 [-0.12, -0.02] |
-0.02 [-0.03, -0.00] |
-0.06 [-0.10, -0.03] |
-0.02 [-0.06, 0.03] |
0.04 [0.03, 0.05] |
0.05 [0.02, 0.08] |
0.05 [0.01, 0.09] |
| Yes |
0.01 [-0.01, 0.03] |
0.03 [-0.04, 0.10] |
-0.04 [-0.13, 0.05] |
-0.11 [-0.14, -0.09] |
-0.09 [-0.16, -0.03] |
-0.12 [-0.21, -0.04] |
-0.02 [-0.04, 0.01] |
-0.01 [-0.07, 0.06] |
0.02 [-0.06, 0.10] |
0.06 [0.04, 0.08] |
0.05 [-0.01, 0.10] |
0.09 [0.02, 0.16] |
|
| School Challenging | -1 SD |
0.05 [0.03, 0.07] |
0.04 [-0.02, 0.11] |
0.05 [-0.03, 0.12] |
-0.10 [-0.12, -0.08] |
-0.13 [-0.19, -0.08] |
-0.04 [-0.11, 0.03] |
-0.03 [-0.05, -0.02] |
-0.07 [-0.13, -0.00] |
-0.06 [-0.12, 0.00] |
0.06 [0.04, 0.07] |
0.11 [0.06, 0.17] |
0.09 [0.04, 0.15] |
| M |
0.06 [0.05, 0.08] |
0.04 [0.00, 0.08] |
0.04 [-0.00, 0.09] |
-0.11 [-0.13, -0.10] |
-0.14 [-0.18, -0.10] |
-0.08 [-0.12, -0.03] |
-0.03 [-0.04, -0.02] |
-0.07 [-0.12, -0.02] |
-0.06 [-0.10, -0.02] |
0.05 [0.04, 0.06] |
0.10 [0.06, 0.14] |
0.05 [0.02, 0.08] |
|
| +1 SD |
0.07 [0.05, 0.09] |
0.04 [-0.02, 0.10] |
0.04 [-0.03, 0.12] |
-0.13 [-0.14, -0.11] |
-0.15 [-0.21, -0.09] |
-0.12 [-0.19, -0.04] |
-0.03 [-0.04, -0.01] |
-0.08 [-0.14, -0.02] |
-0.06 [-0.12, 0.00] |
0.05 [0.04, 0.07] |
0.09 [0.04, 0.14] |
0.01 [-0.05, 0.06] |
|
| School Satisfaction | -1 SD |
0.07 [0.05, 0.09] |
0.09 [0.04, 0.15] |
0.10 [0.02, 0.18] |
-0.10 [-0.12, -0.09] |
-0.01 [-0.07, 0.04] |
-0.07 [-0.14, 0.01] |
-0.03 [-0.05, -0.02] |
-0.01 [-0.06, 0.05] |
-0.11 [-0.18, -0.05] |
0.06 [0.04, 0.07] |
0.05 [0.00, 0.09] |
0.07 [0.01, 0.12] |
| M |
0.06 [0.05, 0.08] |
0.08 [0.04, 0.12] |
0.05 [-0.00, 0.10] |
-0.11 [-0.13, -0.10] |
-0.01 [-0.04, 0.03] |
-0.08 [-0.12, -0.03] |
-0.03 [-0.04, -0.02] |
-0.01 [-0.05, 0.03] |
-0.06 [-0.10, -0.02] |
0.05 [0.04, 0.06] |
0.04 [0.01, 0.08] |
0.05 [0.01, 0.09] |
|
| +1 SD |
0.05 [0.03, 0.07] |
0.06 [0.01, 0.12] |
-0.00 [-0.08, 0.07] |
-0.13 [-0.14, -0.11] |
-0.01 [-0.06, 0.05] |
-0.09 [-0.16, -0.01] |
-0.03 [-0.04, -0.01] |
-0.01 [-0.07, 0.05] |
-0.01 [-0.07, 0.06] |
0.05 [0.04, 0.06] |
0.04 [-0.01, 0.09] |
0.03 [-0.03, 0.09] |
7.4.2 Nonlinear Models
spl_moder_tab_fun <- function(d, weight){
options(papaja.na_string = "")
sup_num <- if(weight == "weighted") 5 else paste0("S",sup_tables)
cap <- sprintf("<strong>Table %s</strong><br><em>Moderators of Pre-Post Change Within-Groups for %s Models", sup_num, weight)
tab <- d %>%
mutate(groups = factor(groups, levels = c("Jail", "CommServ")),
moder = factor(moder, levels = unique(moders2$new_Moder)),
t2 = factor(t2, levels = unique(d$t2))) %>%
arrange(groups, moder, t2) %>%
select(-groups, -t1) %>%
kable(., "html", booktabs = T, escape = F,
col.names = c("Moderator", "Group", rep("b [CI]", times = 4)),
caption = sprintf("<strong>Table X</strong><br><em>", weight),
align = c("l", "l", rep("c", 4))) %>%
kable_styling(full_width = F) %>%
add_header_above(c(" " = 2, "Sensation Seeking" = 1, "Impulsivity" = 1, "Depression" = 1, "Self Esteem" = 1)) %>%
kableExtra::group_rows("Incarceration", 1, 11) %>%
kableExtra::group_rows("Community Service", 12, 22)
save_kable(tab, file = sprintf("%s/results/tables/Table_%s_%s_wi_person_moder_level_tab.html", local_path, sup_num, weight))
if(weight=="weighted") sup_tables <<- sup_tables + 1
return(tab)
}
# tv_group:modvalue, differences in difference in level of contact group after contact as a function of mod
# groups:modvalue, differences in difference in level of contact and no contact groups as a function of mod val
# wave0:tv_group:modvalue, differences in difference in slopes of contact group after contact as a function of mod
# wave0:groups:modvalue, differences in difference in slope of contact and no contact groups as a function of mod val
nl2g_moder_tab <- results.tab %>%
filter(mod != "linear" & grepl(":modval", term) & !grepl("wave", term)) %>%
mutate(term = str_remove_all(term, "groups"),
term = str_remove_all(term, "wave0:"),
term = str_replace_all(term, "modvalue", moder)) %>%
separate(term, c("t1", "t2"), sep = ":") %>%
pivot_wider(names_from = "est", values_from = "value") %>%
mutate(t1 = ifelse(t1 == "tv_group", "Contact: Pre-Post", "None v Contact"),
moder = mapvalues(moder, moders2$Moder, moders2$new_Moder),
t2 = mapvalues(t2, moders3$old, moders3$new),
estimate = sprintf("%s<br>%s", b, CI)) %>%
select(-b, -CI) %>%
pivot_wider(names_from = "pred", values_from = "estimate") %>%
select(groups, weighted, moder, t1, t2, contains(traits$old)) %>%
filter(t1 != "None v Contact") %>%
group_by(weighted) %>%
nest() %>%
ungroup() %>%
mutate(tab = map2(data, weighted, spl_moder_tab_fun))7.4.2.1 Simple Effects
nl2g_moder_sf_tab <- function(d, weight){
options(papaja.na_string = "")
sup_num <- paste0("S",sup_tables)
cap <- sprintf("<strong>Table %s</strong><br><em>Simples Effects of Moderators of Pre-Post Change Within-Groups for %s Models", sup_num, weight)
tab <- d %>%
select(-groups) %>%
kable(.
, "html"
, escape = F
, col.names = c("Moderator", "", rep("b [CI]", 8))
, align = c("r", "r", rep("c", 8))
, caption = cap) %>%
kable_styling(full_width = F) %>%
add_header_above(c(" " = 2, "Pre" = 1, "Post" = 1, "Pre" = 1, "Post" = 1,
"Pre" = 1, "Post" = 1, "Pre" = 1, "Post" = 1)) %>%
add_header_above(c(" " = 2, "Sensation Seeking" = 2, "Impulsivity" = 2,
"Depressive Symptoms" = 2, "Self-Esteem" = 2)) %>%
collapse_rows(1, valign = "top") %>%
kableExtra::group_rows("Community Service", 1, 23) %>%
kableExtra::group_rows("Incarceration", 24, 46)
save_kable(tab, file = sprintf("%s/results/tables/Table_%s_%s_wi_person_moder_simple_effects_tab.html", local_path, sup_num, weight))
sup_tables <<- sup_tables + 1
return(tab)
}
levs <- paste(rep(traits$old, each = 2), c("Pre", "Post"), sep = "_")
nl2g_moder_sf <- nested.moder.mods %>%
filter(mod == "spline") %>%
select(mod, weighted:moder, hypoth) %>%
unnest(hypoth) %>%
separate(term, c("term2", "term", "modval"), sep = "_") %>%
mutate(modval = ifelse(is.na(modval), str_remove_all(term, "[A-Z a-z]"), modval),
modval = mapvalues(modval, c("lower", "upper"), c("-1 SD", "+1 SD")),
modval = paste0(moder, modval),
modval = mapvalues(modval, moders2$comb, moders2$new)) %>%
mutate(Estimate = sprintf("%.2f<br>[%.2f, %.2f]", Estimate, lwr, upr)) %>%
select(weighted:Estimate) %>%
left_join(
nested.moder.mods %>%
filter(mod == "spline") %>%
select(mod, weighted:moder, pool) %>%
unnest(pool) %>%
filter(grepl("tv_group:mod", term)) %>%
mutate(term = str_remove_all(term, "tv_group:"),
sig = ifelse(sign(CI.lower)==sign(CI.upper), "sig", "ns")) %>%
select(weighted:moder, term, sig)
) %>%
mutate(Estimate = ifelse(sig == "sig" & !is.na(sig),
sprintf("<strong>%s</strong>", Estimate), Estimate),
moder = factor(moder, moders2$Moder, moders2$new_Moder)) %>%
select(-sig, -term) %>%
pivot_wider(names_from = c("pred", "term2")
, values_from = "Estimate"
, names_sep = "_") %>%
select(weighted:modval, one_of(levs)) %>%
group_by(weighted) %>%
nest() %>%
ungroup() %>%
mutate(tab = map2(data, weighted, nl2g_moder_sf_tab))| Moderator | b [CI] | b [CI] | b [CI] | b [CI] | b [CI] | b [CI] | b [CI] | b [CI] | |
|---|---|---|---|---|---|---|---|---|---|
| Community Service | |||||||||
| Gender | Male |
2.86 [2.74, 2.99] |
2.64 [2.52, 2.76] |
2.41 [2.29, 2.54] |
2.34 [2.22, 2.46] |
0.76 [0.67, 0.85] |
0.74 [0.66, 0.83] |
3.06 [2.97, 3.15] |
3.14 [3.06, 3.22] |
| Female |
2.53 [2.39, 2.67] |
2.73 [2.59, 2.87] |
2.41 [2.27, 2.55] |
2.30 [2.17, 2.44] |
0.96 [0.85, 1.06] |
0.92 [0.83, 1.02] |
3.10 [3.00, 3.20] |
3.11 [3.01, 3.20] |
|
| Race | White |
2.78 [2.64, 2.93] |
2.77 [2.63, 2.91] |
2.30 [2.15, 2.44] |
2.31 [2.17, 2.45] |
0.84 [0.73, 0.95] |
0.78 [0.67, 0.88] |
3.09 [2.99, 3.19] |
3.09 [3.00, 3.19] |
| Black |
2.50 [2.35, 2.65] |
2.58 [2.43, 2.72] |
2.57 [2.41, 2.72] |
2.38 [2.24, 2.53] |
0.93 [0.80, 1.06] |
0.87 [0.76, 0.98] |
3.15 [3.03, 3.26] |
3.15 [3.05, 3.25] |
|
| Other |
2.83 [2.67, 2.99] |
2.69 [2.54, 2.85] |
2.47 [2.31, 2.63] |
2.28 [2.13, 2.43] |
0.72 [0.59, 0.85] |
0.77 [0.65, 0.89] |
3.01 [2.89, 3.12] |
3.16 [3.05, 3.27] |
|
| Corporal Punishment | No |
2.83 [2.67, 2.98] |
2.74 [2.59, 2.89] |
2.36 [2.21, 2.51] |
2.31 [2.17, 2.45] |
0.83 [0.70, 0.96] |
0.69 [0.57, 0.81] |
3.05 [2.93, 3.18] |
3.22 [3.11, 3.34] |
| Yes |
2.63 [2.49, 2.77] |
2.64 [2.51, 2.77] |
2.47 [2.34, 2.61] |
2.37 [2.24, 2.50] |
0.82 [0.73, 0.91] |
0.85 [0.77, 0.93] |
3.09 [3.01, 3.17] |
3.09 [3.02, 3.17] |
|
| Remedial Education | -1 SD |
2.68 [2.53, 2.82] |
2.69 [2.56, 2.83] |
2.47 [2.33, 2.60] |
2.28 [2.15, 2.41] |
0.76 [0.65, 0.86] |
0.80 [0.70, 0.90] |
3.12 [3.03, 3.22] |
3.23 [3.14, 3.31] |
| M |
2.66 [2.53, 2.79] |
2.65 [2.52, 2.78] |
2.44 [2.31, 2.57] |
2.35 [2.22, 2.47] |
0.83 [0.75, 0.92] |
0.81 [0.73, 0.89] |
3.09 [3.01, 3.17] |
3.14 [3.07, 3.21] |
|
| +1 SD |
2.64 [2.50, 2.79] |
2.61 [2.47, 2.75] |
2.42 [2.28, 2.56] |
2.41 [2.28, 2.55] |
0.91 [0.82, 1.01] |
0.82 [0.73, 0.91] |
3.05 [2.96, 3.13] |
3.06 [2.98, 3.13] |
|
| Parent Died | None |
2.74 [2.59, 2.88] |
2.74 [2.59, 2.88] |
2.42 [2.28, 2.57] |
2.34 [2.19, 2.48] |
0.80 [0.70, 0.91] |
0.79 [0.69, 0.90] |
3.09 [3.00, 3.19] |
3.16 [3.07, 3.26] |
| One |
2.68 [2.53, 2.84] |
2.63 [2.49, 2.77] |
2.41 [2.26, 2.56] |
2.31 [2.18, 2.45] |
0.86 [0.74, 0.97] |
0.81 [0.71, 0.91] |
3.07 [2.96, 3.18] |
3.15 [3.06, 3.24] |
|
| Two |
2.54 [2.21, 2.87] |
2.81 [2.51, 3.10] |
2.52 [2.20, 2.83] |
2.40 [2.12, 2.69] |
1.27 [0.96, 1.58] |
0.84 [0.59, 1.10] |
3.33 [3.05, 3.61] |
2.88 [2.65, 3.10] |
|
| Welfare | No |
2.78 [2.63, 2.92] |
2.74 [2.59, 2.88] |
2.32 [2.18, 2.46] |
2.31 [2.17, 2.46] |
0.81 [0.70, 0.92] |
0.81 [0.70, 0.92] |
3.08 [2.98, 3.18] |
3.18 [3.08, 3.28] |
| Yes |
2.64 [2.51, 2.77] |
2.65 [2.53, 2.78] |
2.51 [2.38, 2.64] |
2.34 [2.22, 2.46] |
0.86 [0.77, 0.96] |
0.81 [0.73, 0.89] |
3.11 [3.02, 3.20] |
3.11 [3.04, 3.19] |
|
| Parental Divorce | No |
2.69 [2.55, 2.82] |
2.62 [2.49, 2.74] |
2.43 [2.30, 2.56] |
2.31 [2.18, 2.43] |
0.84 [0.75, 0.93] |
0.82 [0.74, 0.90] |
3.11 [3.02, 3.19] |
3.15 [3.07, 3.22] |
| Yes |
2.68 [2.54, 2.83] |
2.77 [2.63, 2.90] |
2.38 [2.23, 2.52] |
2.36 [2.22, 2.49] |
0.83 [0.71, 0.94] |
0.79 [0.68, 0.89] |
3.03 [2.93, 3.14] |
3.09 [2.99, 3.18] |
|
| School Challenging | -1 SD |
2.73 [2.58, 2.88] |
2.61 [2.47, 2.75] |
2.34 [2.19, 2.49] |
2.35 [2.21, 2.48] |
0.91 [0.81, 1.01] |
0.87 [0.78, 0.96] |
3.00 [2.90, 3.09] |
3.08 [3.00, 3.17] |
| M |
2.65 [2.51, 2.80] |
2.64 [2.50, 2.78] |
2.35 [2.21, 2.50] |
2.27 [2.13, 2.41] |
0.85 [0.76, 0.94] |
0.82 [0.74, 0.91] |
3.09 [3.01, 3.17] |
3.13 [3.05, 3.21] |
|
| +1 SD |
2.58 [2.41, 2.74] |
2.67 [2.51, 2.83] |
2.37 [2.21, 2.53] |
2.20 [2.04, 2.36] |
0.79 [0.68, 0.90] |
0.78 [0.67, 0.88] |
3.18 [3.08, 3.29] |
3.18 [3.09, 3.28] |
|
| School Satisfaction | -1 SD |
2.57 [2.41, 2.72] |
2.70 [2.54, 2.85] |
2.37 [2.21, 2.52] |
2.23 [2.08, 2.38] |
0.85 [0.74, 0.95] |
0.85 [0.75, 0.96] |
3.13 [3.03, 3.23] |
3.11 [3.01, 3.20] |
| M |
2.65 [2.51, 2.80] |
2.64 [2.50, 2.78] |
2.39 [2.25, 2.53] |
2.30 [2.16, 2.44] |
0.86 [0.77, 0.95] |
0.84 [0.76, 0.92] |
3.08 [2.99, 3.16] |
3.11 [3.04, 3.19] |
|
| +1 SD |
2.74 [2.59, 2.89] |
2.59 [2.45, 2.73] |
2.41 [2.26, 2.55] |
2.37 [2.23, 2.51] |
0.88 [0.78, 0.98] |
0.83 [0.73, 0.92] |
3.02 [2.93, 3.12] |
3.12 [3.04, 3.21] |
|
| Incarceration | |||||||||
| Gender | Male |
2.79 [2.69, 2.89] |
2.67 [2.58, 2.76] |
2.50 [2.40, 2.60] |
2.57 [2.48, 2.66] |
0.71 [0.64, 0.79] |
0.79 [0.72, 0.86] |
3.17 [3.10, 3.24] |
3.17 [3.11, 3.23] |
| Female |
2.63 [2.50, 2.76] |
2.66 [2.53, 2.79] |
2.38 [2.25, 2.50] |
2.61 [2.49, 2.74] |
1.01 [0.90, 1.12] |
1.10 [1.00, 1.20] |
3.07 [2.97, 3.16] |
3.06 [2.98, 3.15] |
|
| Race | White |
2.85 [2.73, 2.97] |
2.77 [2.65, 2.88] |
2.39 [2.27, 2.51] |
2.45 [2.33, 2.56] |
0.85 [0.75, 0.95] |
0.88 [0.79, 0.97] |
3.20 [3.11, 3.28] |
3.14 [3.06, 3.22] |
| Black |
2.62 [2.49, 2.75] |
2.53 [2.41, 2.64] |
2.51 [2.38, 2.64] |
2.64 [2.53, 2.75] |
0.70 [0.58, 0.81] |
0.85 [0.76, 0.95] |
3.09 [3.00, 3.19] |
3.17 [3.09, 3.25] |
|
| Other |
2.79 [2.66, 2.92] |
2.81 [2.68, 2.94] |
2.37 [2.24, 2.50] |
2.64 [2.51, 2.77] |
0.75 [0.64, 0.87] |
0.85 [0.74, 0.96] |
3.13 [3.02, 3.23] |
3.13 [3.04, 3.23] |
|
| Corporal Punishment | No |
2.70 [2.58, 2.82] |
2.67 [2.55, 2.79] |
2.43 [2.31, 2.54] |
2.57 [2.46, 2.68] |
0.78 [0.67, 0.88] |
0.81 [0.72, 0.91] |
3.22 [3.12, 3.31] |
3.16 [3.07, 3.24] |
| Yes |
2.75 [2.64, 2.86] |
2.67 [2.56, 2.77] |
2.48 [2.37, 2.59] |
2.61 [2.51, 2.71] |
0.79 [0.71, 0.87] |
0.88 [0.81, 0.96] |
3.13 [3.06, 3.20] |
3.14 [3.08, 3.21] |
|
| Remedial Education | -1 SD |
2.71 [2.59, 2.82] |
2.63 [2.52, 2.74] |
2.37 [2.26, 2.48] |
2.47 [2.36, 2.57] |
0.79 [0.70, 0.88] |
0.81 [0.72, 0.89] |
3.23 [3.15, 3.31] |
3.22 [3.15, 3.30] |
| M |
2.70 [2.60, 2.80] |
2.63 [2.53, 2.73] |
2.42 [2.32, 2.52] |
2.55 [2.45, 2.64] |
0.77 [0.69, 0.85] |
0.85 [0.78, 0.92] |
3.16 [3.09, 3.23] |
3.16 [3.10, 3.22] |
|
| +1 SD |
2.69 [2.58, 2.80] |
2.64 [2.53, 2.74] |
2.46 [2.36, 2.57] |
2.62 [2.52, 2.72] |
0.75 [0.66, 0.84] |
0.89 [0.81, 0.97] |
3.09 [3.02, 3.17] |
3.10 [3.04, 3.17] |
|
| Parent Died | None |
2.79 [2.67, 2.91] |
2.68 [2.56, 2.81] |
2.41 [2.29, 2.52] |
2.66 [2.54, 2.78] |
0.73 [0.64, 0.82] |
0.89 [0.80, 0.98] |
3.15 [3.06, 3.23] |
3.14 [3.06, 3.21] |
| One |
2.66 [2.54, 2.79] |
2.64 [2.53, 2.75] |
2.54 [2.42, 2.66] |
2.54 [2.44, 2.65] |
0.87 [0.76, 0.97] |
0.85 [0.76, 0.93] |
3.16 [3.06, 3.25] |
3.16 [3.09, 3.24] |
|
| Two |
2.61 [2.31, 2.91] |
2.81 [2.56, 3.07] |
2.84 [2.55, 3.12] |
2.49 [2.25, 2.73] |
1.11 [0.82, 1.41] |
0.87 [0.65, 1.09] |
3.08 [2.85, 3.32] |
3.12 [2.93, 3.30] |
|
| Welfare | No |
2.83 [2.71, 2.95] |
2.77 [2.64, 2.90] |
2.34 [2.22, 2.46] |
2.44 [2.32, 2.56] |
0.80 [0.70, 0.89] |
0.88 [0.78, 0.97] |
3.24 [3.16, 3.33] |
3.19 [3.11, 3.28] |
| Yes |
2.72 [2.61, 2.82] |
2.66 [2.57, 2.75] |
2.46 [2.36, 2.56] |
2.60 [2.51, 2.69] |
0.78 [0.69, 0.86] |
0.86 [0.79, 0.93] |
3.07 [2.99, 3.15] |
3.14 [3.08, 3.20] |
|
| Parental Divorce | No |
2.75 [2.64, 2.85] |
2.60 [2.50, 2.70] |
2.48 [2.38, 2.58] |
2.63 [2.54, 2.72] |
0.84 [0.76, 0.93] |
0.91 [0.83, 0.99] |
3.12 [3.04, 3.19] |
3.09 [3.03, 3.16] |
| Yes |
2.74 [2.63, 2.86] |
2.80 [2.69, 2.91] |
2.39 [2.28, 2.50] |
2.49 [2.39, 2.60] |
0.70 [0.61, 0.80] |
0.81 [0.72, 0.89] |
3.16 [3.07, 3.24] |
3.21 [3.14, 3.28] |
|
| School Challenging | -1 SD |
2.67 [2.55, 2.79] |
2.62 [2.51, 2.73] |
2.43 [2.32, 2.55] |
2.56 [2.46, 2.67] |
0.82 [0.72, 0.91] |
0.86 [0.78, 0.94] |
3.10 [3.02, 3.18] |
3.14 [3.07, 3.21] |
| M |
2.67 [2.55, 2.78] |
2.60 [2.49, 2.70] |
2.45 [2.34, 2.55] |
2.58 [2.48, 2.68] |
0.76 [0.68, 0.84] |
0.82 [0.75, 0.90] |
3.15 [3.08, 3.22] |
3.16 [3.09, 3.22] |
|
| +1 SD |
2.66 [2.54, 2.79] |
2.57 [2.45, 2.69] |
2.46 [2.34, 2.58] |
2.59 [2.47, 2.71] |
0.71 [0.61, 0.80] |
0.78 [0.69, 0.88] |
3.20 [3.11, 3.28] |
3.17 [3.09, 3.25] |
|
| School Satisfaction | -1 SD |
2.66 [2.54, 2.78] |
2.68 [2.57, 2.80] |
2.44 [2.32, 2.56] |
2.55 [2.44, 2.66] |
0.73 [0.64, 0.83] |
0.86 [0.77, 0.94] |
3.13 [3.05, 3.22] |
3.16 [3.08, 3.23] |
| M |
2.68 [2.57, 2.79] |
2.62 [2.51, 2.72] |
2.44 [2.33, 2.55] |
2.57 [2.47, 2.67] |
0.77 [0.69, 0.84] |
0.84 [0.76, 0.91] |
3.15 [3.08, 3.22] |
3.15 [3.09, 3.22] |
|
| +1 SD |
2.69 [2.58, 2.81] |
2.55 [2.44, 2.66] |
2.44 [2.33, 2.55] |
2.58 [2.48, 2.69] |
0.80 [0.71, 0.89] |
0.81 [0.74, 0.89] |
3.16 [3.08, 3.24] |
3.15 [3.08, 3.22] |
|
7.5 Plot Results
7.5.1 Linear
7.5.2 Non-Linear
plot_moder_fun <- function(d, moder, weight, mod, gr){
sup_num <- paste0("S",sup_figs)
md <- moders2 %>% filter(Moder == moder)
cols <- c("black", "grey40", "grey70")[1:nrow(md)]
if(mod == "spline"){
p <- d %>%
mutate(wave0 = ifelse(abs(wave0) > .05, round(wave0, 1), wave0)) %>%
distinct() %>%
full_join(d %>% select(dv, groups, modvalue) %>%
filter(groups != "None") %>%
distinct() %>%
mutate(wave0 = 0, pred = NA)) %>%
select(-tv_group, -age_c) %>%
group_by(dv, groups, modvalue, wave0) %>%
summarise_all(mean, na.rm = T) %>%
ungroup() %>%
distinct() %>%
mutate(modvalue = mapvalues(modvalue, md$old, md$new)
, dv = factor(dv, traits$old, traits$breaks)) %>%
ggplot(aes(x = wave0
, y = pred
, color = modvalue
, linetype = modvalue
, group = modvalue)) +
scale_color_manual(values = cols) +
scale_x_continuous(limits = c(-10.2, 10.2), breaks = seq(-10, 10, 5)) +
geom_line(size = 1) +
labs(x = "Wave (Centered at Average Age of Contact)"
, y = "Predicted Level"
, color = NULL
, linetype = NULL
, title = sprintf("Figure %s. %s Nonlinear Trajectories Across Contact Groups and %s"
, sup_num, str_to_title(weight), unique(md$new_Moder))) +
facet_grid(groups ~ dv) +
theme_classic() +
theme(panel.background = element_rect(color = "black", size = .8)
, strip.background = element_blank()
, strip.text = element_text(face = "bold", size = rel(1.1))
, legend.position = "bottom"
, legend.text = element_text(face = "bold", size = rel(.9))
, axis.line = element_line(size = .1)
, axis.text = element_text(face = "bold", size = rel(1), color = "black")
, axis.title = element_text(face = "bold", size = rel(1.1)),
plot.title = element_text(face = "bold", size = rel(1.2), hjust = .5))
} else {
p <- d %>%
group_by(dv, age0, groups, modvalue) %>%
summarize_all(mean) %>%
mutate(dv = factor(dv, levels = traits$old, traits$breaks),
groups = factor(groups, levels = c("None", "Jail", "CommServ")
, labels = c("Neither", "Incarceration", "Community\nService"))
, modvalue = mapvalues(modvalue, md$old, md$new)) %>%
ggplot(aes(x = age0*5 + 14
, y = pred
, color = modvalue
, linetype = modvalue
, group = modvalue)) +
scale_color_manual(values = cols) + # colors
geom_line(size = .9) + # add the line
labs(x = "Age (Centered at 14)" # titles
, y = "Predicted Level"
, color = NULL
, linetype = NULL
, title = sprintf("Figure %s. %s Linear Trajectories Across Contact Groups and %s"
, sup_num, str_to_title(weight), unique(md$new_Moder))) +
facet_grid(groups ~ dv) + # faceting
theme_classic() +
theme(panel.background = element_rect(color = "black", size = .8),
axis.text = element_text(face = "bold", size = rel(1.2)),
axis.title = element_text(face = "bold", size = rel(1.2)),
strip.text = element_text(face = "bold", size = rel(1.2)),
strip.background = element_blank(),
legend.text = element_text(face = "bold"),
legend.title = element_text(face = "bold", size = rel(1.2)),
legend.position = "bottom",
plot.title = element_text(face = "bold", size = rel(1.2), hjust = .5))
}
ggsave(p
, file = sprintf("%s/results/plots/Figure_%s_%s_%s_trajectories_%s.png", local_path, sup_num, mod, weight, moder)
, height = 6
, width = 10)
sup_figs <<- sup_figs + 1
return(p)
}
moder_plots <- nested.moder.mods %>%
group_by(mod) %>%
nest() %>%
ungroup() %>%
mutate(data = map(data, ~(.) %>%
select(weighted, dv = pred, mod_type = groups, moder, pred.fx) %>%
unnest(pred.fx) %>%
group_by(weighted, moder, mod_type) %>%
nest() %>%
ungroup())) %>%
unnest(data) %>%
mutate(pmap(list(data, moder, weighted, mod, mod_type), plot_moder_fun))
6 Socialization Effects: Growth Models
6.1 Zero Order Relationships
First, we look at the correlations between each of the 4 outcomes across ages.
6.2 Measurement Invariance
6.3 Basic Models
Next, we fit the growth models with and without weights for each of the 4 outcomes (CESD, Impulsivity, Sensation Seeking, and Self Esteem) using a linear of the following forms.
<strong?Linear:Then, to see if the timing of the event has a substantive impact on trajectories of change, we will run a series of not-preregistered non-linear multilevel spline models.
First, to compare whether the trajectory of those who experience either incarceration and community service differs from those who don’t (i.e. a change point), we will run the following model:-
Level 1: \(Y_{ij} = \beta_{0j} + \beta_{1j}*time0_{ij} + \beta_{2j}*D_{ij} + \beta_{3j}*time0_{ij}*D_{ij} + \varepsilon{ij}\)
-
Level 2:
-
\(\beta_{0j} = \gamma_{00} + \gamma_{01}*age_j + U_{0j}\)
-
\(\beta_{1j} = \gamma_{10} + + \gamma_{01}*age_j + U_{1j}\)
-
\(\beta_{2j} = \gamma_{20}\)
-
\(\beta_{3j} = \gamma_{30}\)
Second, to compare (1) whether any contact with the criminal justice system impacts trajectories of change and (2) whether trajectories differ depending on court ordered community service versus incarceration, we will the following model:6.3.1 Define Groups
6.4 Setup Data
6.4.1 Basic Linear Growth Model
6.4.2 Two-Group Model
6.5 Define Functions
6.5.1 Model Functions
6.5.2 Pool Models
Next, we pool results across models. Below is a longish function for pooling the fixed and random effects of the models of the imputed datasets.
6.5.3 Pooled Random Effects
For plotting individual trajectories and effect sizes, we need random effects. So we pool those across models
6.5.4 Pooled
vcov()MatrixTo plot confidence bands around trajectories, we’re going to need the variance-covariance matrix. So we pool those across imputations.
6.5.5 Standard Error Function
To unpack the cross-level interactions, we need the standard errors. We use the method from Preacher and Curran (2006).
6.5.6 Hypothesis Tests
We want to get individual slope estimates for each Group for each trait. To get the confidence intervals around them in addition to their estimates, we use the
glht()function from themultcomppackage.6.6 Run Models
6.7 Load Results
Now, we actually pool the fixed and random effects, both for average effects, as well for individual participants, using the functions we defined above.
6.8 Tables
6.8.1 Linear Slopes
6.8.2 Nonlinear Levels
Weighted Estimated Levels Across Groups at Average Contact Points and Differences Between Groups and Across Contact Points
b [CI]
b [CI]
b [CI]
b [CI]
b [CI]
b [CI]
b [CI]
b [CI]
[2.51, 2.55]
[0.06, 0.22]
[2.20, 2.24]
[0.28, 0.44]
[0.65, 0.68]
[0.14, 0.25]
[3.25, 3.27]
[-0.16, -0.06]
[2.66, 2.82]
[-0.10, -0.03]
[2.37, 2.52]
[0.10, 0.17]
[0.72, 0.84]
[0.05, 0.11]
[3.10, 3.20]
[-0.03, 0.03]
[2.59, 2.75]
[2.51, 2.65]
[0.81, 0.91]
[3.10, 3.20]
[2.53, 2.57]
[0.03, 0.23]
[2.18, 2.22]
[0.03, 0.23]
[0.64, 0.66]
[0.10, 0.22]
[3.25, 3.28]
[-0.19, -0.08]
[2.59, 2.80]
[-0.06, 0.02]
[2.32, 2.52]
[-0.13, -0.05]
[0.77, 0.90]
[-0.06, 0.01]
[3.02, 3.15]
[0.01, 0.08]
[2.58, 2.78]
[2.23, 2.43]
[0.75, 0.87]
[3.07, 3.18]
6.8.3 Nonlinear Slopes
Weighted Estimated Levels Across Groups at Average Contact Points and Differences Between Groups and Across Contact Points
b [CI]
b [CI]
b [CI]
b [CI]
b [CI]
b [CI]
b [CI]
b [CI]
[0.01, 0.01]
[-0.02, -0.01]
[-0.02, -0.01]
[-0.03, -0.01]
[-0.01, -0.01]
[-0.02, -0.01]
[0.00, 0.01]
[-0.00, 0.01]
[0.02, 0.04]
[-0.04, -0.02]
[-0.02, -0.01]
[-0.03, -0.02]
[-0.02, -0.00]
[-0.01, -0.00]
[0.01, 0.02]
[-0.01, 0.00]
[-0.01, 0.01]
[-0.04, -0.03]
[-0.02, -0.01]
[0.00, 0.01]
[0.01, 0.01]
[-0.01, 0.00]
[-0.02, -0.01]
[-0.01, 0.01]
[-0.01, -0.01]
[0.00, 0.01]
[0.00, 0.01]
[0.00, 0.01]
[0.02, 0.04]
[-0.03, -0.02]
[-0.01, 0.01]
[-0.02, -0.01]
[-0.02, -0.00]
[0.01, 0.02]
[-0.01, 0.01]
[0.00, 0.02]
[-0.00, 0.01]
[-0.02, -0.01]
[-0.00, 0.01]
[0.01, 0.02]
6.9 Results Plots
6.9.1 Fixed Effects
Next, we plot the results by extracting the coefficients from the models and using the simple slopes for each of the groups.
6.9.2 Trajectories
6.9.2.1 Linear
For example, the simple slopes for the fixed effects linear model can be defined as follows (Note the absence of the random terms in the simple slopes).
We extract each of the coefficients from the pooled models and multiply them by the expanded model frame (all combinations of plausible levels of groups and age).
6.9.2.2 Non-Linear
6.9.3 Random Effects